To find the intersection of two great circles defined by the arcs
from pt1 (lat1, lon1) to pt2 and from pt3 to pt4.
-----------------
Use the "W longitude is positive" convention of the Aviation
Formulary. (For the more conventional convention, change the sign of
all longitudes in the following.)
For each point we can associate a unit vector pointing to it from
the center of the earth whose components are:
e ={ex,ey,ez} = {cos(lat)*cos(lon), -cos(lat)*sin(lon), sin(lat)} (1)
which we can invert with
lat=atan2(ez, sqrt(ex^2 + ey^2)); lon=atan2(-ey, ex) (2)
For any great circle, defined by pts 1 and 2, the point
P(e1,e2) =(e1 X e2)/||e1 X e2||
is perpendicular the the plane of the circle. (Its negative is the opposite
point).
Here e1 X e2 is the vector cross-product whose components are
{e1y *e2z -e2y *e1z, e1z *e2x -e2z *e1x, e1x *e2y -e1y *e2x} (3)
respectively (see later for a numerically robust way to compute this!)
||e|| is the length of a vector defined by
||e|| =sqrt(ex^2 + ey^2 + ez^2) (4)
The intersections of the great circles can be seen to be given by
+- P( P(e1,e2), P(e3,e4) )
Direct computation of the cross-product will fail at small distances
because of rounding error. Application of some trig identities
gives:
e1 X e2 = { sin(lat1-lat2) *sin((lon1+lon2)/2) *cos((lon1-lon2)/2) -
sin(lat1+lat2) *cos((lon1+lon2)/2) *sin((lon1-lon2)/2) ,
sin(lat1-lat2) *cos((lon1+lon2)/2) *cos((lon1-lon2)/2) +
sin(lat1+lat2) *sin((lon1+lon2)/2) *sin((lon1-lon2)/2) ,
cos(lat1)*cos(lat2)*sin(lon1-lon2) } (5)
which avoids this problem.
Algorithm:
compute e1 X e2 and e3 X e4 using (5).
Normalize ea= (e1 X e2)/ ||e1 X e2|| , eb=(e3 x e4)/||e3 X e4|| using (4)
Compute ea X eb using (3)
Invert using (2) (it's unnecessary to normalize first).
The two candidate intersections are (lat,lon) and the antipodal point
(-lat, lon+pi)