Next order the two points such that p1 <= p2.
Next put the two points in cartesial coordinates (x,y,z).
Through a little geometry it can be found that:
(x1,y1,z1) =
(cos(p1)*sin(q1), sin(p1)*sin(q1), cos(q1)) and
(x2,y2,z2) =
(cos(p2)*sin(q2), sin(p2)*sin(q2), cos(q2)).
The angle between two vectors is defined as cos(t)=(v dot w)/(|v| * |w|).
Let v donate the vector from (0,0,0) to (x1,y1,z1) and
let w donate the vector from (0,0,0) to (x2,y2,z2).
The dot product of v and w is:
cos(p1)*sin(q1)*cos(p2)*sin(q2) +
sin(p1)*sin(q1)*sin(p2)*sin(q2) +
cos(q1)*cos(q2).
If we rotate the longitude of both points by p1 then p1 becomes 0 and
p2 becomes p2-p1, in which case the dot product is:
cos(0)*sin(q1)*cos(p2-p1)*sin(q2) +
sin(0)*sin(q1)*sin(p2-p1)*sin(q2) +
cos(q1)*cos(q2). =
sin(q1)*cos(p2-p1)*sin(q2) +
cos(q1)*cos(q2).
Now, both |u| and |v| must be 1 since they lie on the surface of a unit sphere. Thus the angle between the two points must be cos-1[sin(q1)*cos(p2-p1)*sin(q2) + cos(q1)*cos(q2) / (1 * 1). = sin(q1)*cos(p2-p1)*sin(q2) + cos(q1)*cos(q2)].
This angle is the ratio of the distance along the surface between the
two points to the circumference of the earth. Thus the distance is the angle times
the circumference of the earth which equals:
cos-1[(sin(q1)*cos(p2-p1)*sin(q2) +
cos(q1)*cos(q2))] * r.
Michael Shackleford, A.S.A.