|
Also known as: Stupid Geometry Which I Tend To Forget.
The problem is a follows: find the length of the shortest path connecting
two points on a sphere.
Suppose the location of these points are given as latitude/longitude pairs
(b,l).
Step 1: transform those coordinates into Cartesian coordinates.
(x,y,z) = r ( cos b cos l, cos b sin l, sin b)
Step 2: make use of the inner product rule: r1 * r2 = |r1||r2| cos (p)
(where p is the angle between the two vectors r1 and r2)
( cos b1 cos l1, cos b1 sin l1, sin b1) * ( cos b2 cos l2, cos b2 sin l2, sin b2) = cos a
cos b1 cos b2 cos l1 cos l2 + cos b1 cos b2 sin l1 sin l2 + sin b1 sin b2 = cos a
cos b1 cos b2 ( sin l1 sin l2 + cos l1 cos l2) + sin b1 sin b2 = cos a
Step 3: Apply brains. (Well, dig up some oldskool
books containing trigonometric
formulas and do some bookkeeping)
cos (l1-l2) * sin b1*sin b2 + cos b1 * cos b2 = cos a
This is the result that MathWorld gives.
However, cosines are teh suck, and the computational cost can be reduced
even further if it is taken into account that products of sines and cosines
can be written as:
sin a * sin b = ((sin
a-b)-(sin a+b))/2 and cos a*cos b=((sin a-b)+(sin a+b))/2
Now, the thing can be reduced to three cosine lookups. Yay!
cos (l1-l2) * (sin (b1-b2)-(sin b1+b2) + (sin (b1-b2)+sin (b1+b2)) = 2 cos a
We're almost there: the only thing left to calculate the arc length is
to multiply the angle we just found (in radians!) by the Earth's radius,
which is about 6371 km. (No need to be too accurate with the radius,
the Earth is no perfect sphere either.)
Testcase: from my house (52.0103N, 4.3661E) to the faculty of Aerospace
Engineering (51.9897N, 4.3759E) gives 2.39 km as a result. The software in
my GPS handheld agrees with that distance, so I'm happy.
And the moment you've all been waiting for, an implementation in python:
#!/usr/bin/env python
from math import cos,acos,radians,degrees
def lbtoangle(b1,l1,b2,l2):
p1 = cos(radians(l1-l2))
p2 = cos(radians(b1-b2))
p3 = cos(radians(b1+b2))
return degrees(acos(((p1*(p2+p3))+(p2-p3))/2))
print radians(lbtoangle(52.0103,4.3661,51.9897,4.3759))*6371
|