The solution for an ellipsoid is pretty messy--it is an irregular shape, not a circle--and is best computed numerically rather than with a formula.
On a world map the difference between the WGS84 solution and a purely spherical solution will only barely be noticeable (it's about one pixel on a screen). The same difference would be created by changing the mask angle by about 0.2 degrees or in using a polygonal approximation. If these errors are acceptably small, then you can exploit the symmetry of the sphere to obtain a simple formula.
This map (using an Equirectangular projection) shows the coverage for a satellite at 22,164 kilometers (from the earth's center) with a mask angle of m = 15 degrees on the WGS84 spheroid. Recomputing the coverage for a sphere does not visibly change this map.
On the sphere, the coverage will truly be a circle centered at the satellite's location, so we only need to figure out its radius, which is an angle. Call this t. In cross-section there is a triangle OSP formed by the earth's center (O), the satellite (S), and any point (P) on the circle:
The side OP is the radius of the earth, R.
The side OS is the height of the satellite (above the earth's center). Call this h.
The angle OPS is 90 + m.
The angle SOP is t, which we want to find.
Because the three angles of a triangle sum to 180 degrees, the third angle OSP must equal 90 - (m + t).
The solution is now a matter of elementary trigonometry. The (planar) law of sines asserts that
sin(90 - (m+t)) / r = sin(90 + m) / h.
The solution is
t = ArcCos(cos(m) / (h/r)) - m.
As a check, consider some extreme cases:
When m = 0, t = ArcCos(r/h), which can be verified with elementary Euclidean geometry.
When h = r (the satellite hasn't launched), t = ArcCos(cos(m) / 1) - m = m - m = 0.
When m = 90 degrees, t = ArcCos(0) - 90 = 90 - 90 = 0, as it should be.
This reduces the problem to drawing a circle on the sphere, which can be solved in many ways. For instance, you can buffer the satellite's location by t * R * pi/180 using an equidistant projection centered at the satellite. Techniques for working with circles on the sphere directly are illustrated at https://gis.stackexchange.com/a/53323/664.
Edit
FWIW, for GPS satellites and small mask angles (less than 20 degrees or so), this non-trigonometric approximation is accurate (to a few tenths of a degree and less than a few hundredths of a degree when the mask angle is under 10 degrees):
t (degrees) = -0.0000152198628163333 * (-5.93410042925107*10^6 +
3.88800000000000*10^6 r/h + 65703.6145507725 m +
9.86960440108936 m^2 - 631.654681669719 r/h m^2)
For example, with a mask angle of m = 10 degrees and a satellite at 26,559.7 km above the earth's center (which is the nominal distance of a GPS satellite), this approximation gives 66.32159..., whereas the value (correct for the sphere) is 66.32023... .
(The approximation is based on a Taylor series expansion around m = 0, r/h = 1/4.)
Best Answer
Assuming that you do not worry about the errors produced by treating geographic coordinates as spherical coordinates instead of ellipsoidal ones.
When applying these formulas you are obtaining two planes, one plane for the latitudes and another plane for the longitudes. And for each of them, a couple of Cartesian coordinates.
In the plane of the latitudes, all of them are being treated as if they belonged to the same meridian. In the plane of the longitudes, all of them are being treated as if they belonged to the same parallel.
This is because you are separating the spherical coordinates into two two-dimensional spaces with polar coordinates. Where the radius is the common coordinate to both. That is, from three spherical coordinates [(radius, latitude, longitude)], you are obtaining two pairs of polar coordinates [(radius, latitude) ; (radius, longitude)], and from there converting to two pairs of Cartesian coordinates [(x, y) ; (x, y)].
Now suppose you have a point located at
longitude = 0andlatitude = 0. And another point located in the North pole (assumelongitude = 0andlatitude = 90). And you want to know the position of the 1/3 of the line between them. The distance between parallels is the same in the sphere, regardless of the longitude they are. Therefore, the point will be inlongitude = 0andlatitude = 30degrees.But, what happens if you interpolate them in a Cartesian system (transforming only the latitudes from a sphere with
radius = 1). One point is inx = 1,y = 0. The other is inx = 0,y = 1. The interpolated point is inx = 2/3,y = 1/3. And the arctangent of((1/3) / (2/3))is not 30 degrees. Draw them and you will understand the problem.That method can work for the aritmethic mean between two positions. But as soon as you try to average more than two positions, the result becomes wrong again.
From https://en.wikipedia.org/wiki/Mean_of_circular_quantities#Mean_of_angles:
For the longitudes, the calculation proposed will be worst. Because the distance between meridians is a function of the latitude in wich the distance is measured. Two points separated by one degree of longitude are separated by a different distance according to what latitude they are. That is, longitudes can not be treated as if they all belong to the same latitude.
Therefore, the correct way to handle this depends a little on the particular conditions and needs. As a general rule, the conversion to the geocentric system of Cartesian coordinates works well for locations close to each other. But in the same 3D system, not in two 2D systems.
For means of more than two locations, or interpolations, of points that are not relatively close to each other, we will have to sit down with some books in hand to analyze the best solution for the particular need that arises.