Drawing ground station coverage circles

Question

I'm trying to make a little application to draw ground station coverage circles (before anyone suggests using other software, I enjoy making these kinds of things myself!)

I already have code that converts the lat and long coordinates of a point on the surface to points on a Mercator projection map so now I'm drawing the 'circles' around these.

I also have the equation to calculate the slant range (maximum distance between the satellite and the ground station at minimum elevation).

Here's what I do currently:

1. Convert lat, long and Re of ground stations (spherical coords) to x1, y1, z1 (Cartesian coords)

2. Convert elevation, azimuth and slant range (spherical coords) to x2, y2, z2 (Cartesian coords) in ground station centred reference frame

3. Translate x2,y2,z2 from part 2 to Earth centred (x2+x1, y2+y1, z2+z1) lets call these x3, y3, z3

4. Rotate x3, y3, z3 to align with earth centred reference frame lets say x4, y4, z4

5. Convert x4, y4, z4 to spherical coordinates and then plot the lat and long that I've found.

To my eye the above should give me the correct points to plot however there is an error somewhere. I know I step 1 is correct (the magnitude of the vector is Re). I know step 2 is correct since the magnitude of the vector is the slant range. But here's where I hit a problem.... No idea where the problem is though.

Answer 1

Here's how I would solve it:

First, the math is easier if proper spherical coordinates $(r, \theta, \phi)$ are used instead of longitude and latitude are used. The difference is latitude is the angle from the equator, whereas $\phi$ is the angle from the north pole. To find $\theta$ and $\phi$ from the longitude and latitude, $$\theta = \text{Long}$$ $$\phi = 90 - \text{Lat}$$ All angles are in degrees.

1. Convert the elevation and azimuth of the satellite into absolute angles $\theta$ and $\phi$. This will depend on both the position of the satellite in the sky, and the location of ground station. $$\phi_\text{satellite} = \phi_\text{station} - (90-\text{Elev}) \cos(\frac{2 \pi}{360} \text{Azim})$$ $$\theta_\text{satellite} = \theta_\text{station} + (90-\text{Elev}) \sin(\frac{2 \pi}{360} \text{Azim})$$
2. Convert the vector between the center of the earth and the ground station, and the vector between the station and the satellite to Cartesian coordinates $$x = \text{Range} \cos(\frac{2 \pi}{360}\theta) \sin(\frac{2 \pi}{360}\phi)$$ $$y = \text{Range} \sin(\frac{2 \pi}{360}\theta) \sin(\frac{2 \pi}{360}\phi)$$ $$z = \text{Range} \cos(\frac{2 \pi}{360}\phi)$$ For vector from the center of the earth to the ground station, $\text{Range}$ is the radius of the earth.
3. Add the two vectors together, convert back to spherical coordinates, and plot them. $$r = \sqrt{x^2 + y^2 + z^2}$$ $$\theta = \frac{360}{2 \pi}\tan^{-1}(y/x)$$ $$\phi = \frac{360}{2 \pi}\cos^{-1}(z/r)$$

Enjoy!