Questions like this can often be answered by solving triangles using laws of spherical trigonometry. You can look these up: they are referenced below. Normally, when you know three of the six parts of a spherical triangle, the laws let you find the other three in terms of the sines and cosines of the parts you know. The trick is to draw a useful triangle. When you're dealing with meridians, think of using a pole (which is common to all meridians) as one of the triangle's vertices. This suggests focusing on the triangles N-M-T1 or, equivalently, N-M-T2. Here are the details, a general solution, and a worked example.
Let the point in question be at latitude phi and the circle's radius be r (expressed as an angle, not as a distance. The angle in radians is the radial distance divided by the sphere's radius). We start by finding the difference between the local longitude and the extreme longitudes; let this difference be lambda.
Because the extreme meridians are tangent to the circle, angles M-T1-N and M-T2-N are right angles. This is nice because the sine and cosine of a right angle are easy to find and simple (they equal 1 and 0, respectively). The spherical Law of Sines says
sin(lambda) / sin(r) = sin(90 degrees) / sin(90 - phi)
= 1 / cos(phi).
Solving this gives
lambda = ArcSin( sin(r) / cos(phi) ).
For example, let r = 45 degrees and phi = -30 degrees. The solution is lambda = 54.7356 degrees.
This much you already found out: it's the formula quoted in the question. Now we can apply the spherical Law of Cosines to relate the latitude phi' of the extremal points T1 and T2 to the radius r: this is what we are looking for.
sin(phi) = sin(phi') cos(r).
The solution for the common latitude of the tangent points is
phi' = ArcSin( sin(phi) / cos(r) )
(for a circle centered at latitude phi of radius r radians.)
The remainder of the reply discusses and illustrates this result.
Notice in particular that when we start at the equator, phi = 0, easily giving phi' = 0 as expected. Otherwise, the sine of phi' is larger (in size) than the sine of phi, implying the solution is closer to the nearest pole, again as expected. In the example, phi' works out to -45 degrees.
The south pole is near the bottom of this pseudo-3D image of a sphere. The meridians start at the circle's center and go in 15 degree increments to either side, showing +-15, +-30, and +-45 degree longitudinal displacements, then stop at the limiting values of +-54.7356 degrees. The red dots are situated along these limiting meridians at latitudes of -45 degrees.
These formulas work for any circle that does not include either pole.
In this image, the circle's center (black dot) is just north of 60 degrees north latitude. Approximate circles of latitude are shown as dashed curves in 10 degree increments. The radius, equal almost to 30 degrees of arc, takes this circle almost to the north pole (where the meridians are converging at the top). The extremal meridians are therefore spread far apart (lambda is about 70 degrees) and, accordingly, the points of tangency are also close to the north pole, near 80 degrees latitude. This shows why the points of tangency (red dots) are usually closer to the nearest pole than the circle's center is.
In the Specification for the Wide Area Augmentation System (WAAS) on page 99 there are two formulae for calculating the "pierce point" of a GPS signal through the ionosphere. This may not be your application area, but may be able to adapt these to your application and they happen to work with longitude, latitude, azimuth, and elevation angle. Here's a picture from the document:
On the site http://paulbourke.net/geometry/circlesphere/ there is a section called "Intersection of a Line and a Sphere (or circle)" with more general equations for intersections. At the end of that section you will find some implementations in C, VB, Python, and LISP.
As a bonus (this is a limited time offer, you are one of a select few to receive this offer ;-) ) here is a short Python script for visualizing line-sphere intersections in Blender (adapted from this blog post. Note that is takes Earth-centered, Earth-fixed coordinates, not longitude/latitude/azimuth/elevation angles.
import bpy
EARTH_EQUATOR = 6378137.0 # meters
IONO_HEIGHT = 350000.0 # meters
# scale down the radius by 1E6 so that it fits in blender's world
radius = (EARTH_EQUATOR + IONO_HEIGHT) / 1E6
bpy.ops.mesh.primitive_uv_sphere_add(size=radius)
from mathutils import Vector
# weight
w = 10
def MakePolyLine(objname, curvename, cList):
curvedata = bpy.data.curves.new(name=curvename, type='CURVE')
curvedata.dimensions = '3D'
objectdata = bpy.data.objects.new(objname, curvedata)
objectdata.location = (0,0,0) #object origin
bpy.context.scene.objects.link(objectdata)
polyline = curvedata.splines.new('POLY')
polyline.points.add(len(cList)-1)
for num in range(len(cList)):
polyline.points[num].co = (cList[num])+(w,)
polyline.order_u = len(polyline.points)-1
polyline.use_endpoint_u = True
# here's a sample vector
listOfVectors = [(0.81,3.443,-5.723),(3.723,0.111,-5.603)]
MakePolyLine("NameOfMyCurveObject", "NameOfMyCurve", listOfVectors)
With more vectors, this would output something like the image below:
Best Answer
It's not much harder on the sphere than on the plane, once you recognize that
The points in question are the mutual intersections of three spheres: a sphere centered beneath location x1 (on the earth's surface) of a given radius, a sphere centered beneath location x2 (on the earth's surface) of a given radius, and the earth itself, which is a sphere centered at O = (0,0,0) of a given radius.
The intersection of each of the first two spheres with the earth's surface is a circle, which defines two planes. The mutual intersections of all three spheres therefore lies on the intersection of those two planes: a line.
Consequently, the problem is reduced to intersecting a line with a sphere, which is easy.
Here are the details. The inputs are points P1 = (lat1, lon1) and P2 = (lat2, lon2) on the earth's surface, considered as a sphere, and two corresponding radii r1 and r2.
Convert (lat, lon) to (x,y,z) geocentric coordinates. As usual, because we may choose units of measurement in which the earth has a unit radius,
In the example, P1 = (-90.234036 Degree, 37.673442 Degree) has geocentric coordinates x1 = (-0.00323306, -0.7915, 0.61116) and P2 = (-90.953669 Degree, 36.109997 Degree) has geocentric coordinates x2 = (-0.0134464, -0.807775, 0.589337).
Convert the radii r1 and r2 (which are measured along the sphere) to angles along the sphere. By definition, one nautical mile (NM) is 1/60 degree of arc (which is pi/180 * 1/60 = 0.0002908888 radians). Therefore, as angles,
The geodesic circle of radius r1 around x1 is the intersection of the earth's surface with a Euclidean sphere of radius sin(r1) centered at cos(r1)*x1.
The plane determined by the intersection of the sphere of radius sin(r1) around cos(r1)*x1 and the earth's surface is perpendicular to x1 and passes through the point cos(r1)x1, whence its equation is x.x1 = cos(r1) (the "." represents the usual dot product); likewise for the other plane. There will be a unique point x0 on the intersection of those two planes that is a linear combination of x1 and x2. Writing x0 = ax1 + b*x2 the two planar equations are
Using the fact that x2.x1 = x1.x2, which I shall write as q, the solution (if it exists) is given by
In the running example, I compute a = 0.973503 and b = 0.0260194.
Evidently we need q^2 != 1. This means that x1 and x2 can be neither the same point nor antipodal points.
Now all other points on the line of intersection of the two planes differ from x0 by some multiple of a vector n which is mutually perpendicular to both planes. The cross product
does the job provided n is nonzero: once again, this means that x1 and x2 are neither coincident nor diametrically opposite. (We need to take care to compute the cross product with high precision, because it involves subtractions with a lot of cancellation when x1 and x2 are close to each other.) In the example, n = (0.0272194, -0.00631254, -0.00803124).
Therefore, we seek up to two points of the form x0 + t*n which lie on the earth's surface: that is, their length equals 1. Equivalently, their squared length is 1:
The term with x0.n disappears because x0 (being a linear combination of x1 and x2) is perpendicular to n. The two solutions easily are
and its negative. Once again high precision is called for, because when x1 and x2 are close, x0.x0 is very close to 1, leading to some loss of floating point precision. In the example, t = 1.07509 or t = -1.07509. The two points of intersection therefore equal
Finally, we may convert these solutions back to (lat, lon) by converting geocentric (x,y,z) to geographic coordinates:
For the longitude, use the generalized arctangent returning values in the range -180 to 180 degrees (in computing applications, this function takes both x and y as arguments rather than just the ratio y/x; it is sometimes called "ATan2").
I obtain the two solutions (-88.151426, 36.989311) and (-92.390485, 38.238380), shown in the figure as yellow dots.
The axes display the geocentric (x,y,z) coordinates. The gray patch is the portion of the earth's surface from -95 to -87 degrees longitude, 33 to 40 degrees latitude (marked off with a one degree graticule). The earth's surface has been made partly transparent to show all three spheres. The correctness of the computed solutions is evident by how the yellow points sit at the intersections of the spheres.