Next order the two points such that p_{1} <= p_{2}.

Next put the two points in cartesial coordinates (x,y,z).

Through a little geometry it can be found that:

(x_{1},y_{1},z_{1}) =
(cos(p_{1})*sin(q_{1}), sin(p_{1})*sin(q_{1}), cos(q_{1})) and

(x_{2},y_{2},z_{2}) =
(cos(p_{2})*sin(q_{2}), sin(p_{2})*sin(q_{2}), cos(q_{2})).

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 (x_{1},y_{1},z_{1}) and

let w donate the vector from (0,0,0) to (x_{2},y_{2},z_{2}).

The dot product of v and w is:

cos(p_{1})*sin(q_{1})*cos(p_{2})*sin(q_{2}) +
sin(p_{1})*sin(q_{1})*sin(p_{2})*sin(q_{2}) +
cos(q_{1})*cos(q_{2}).

If we rotate the longitude of both points by p_{1} then p_{1} becomes 0 and
p_{2} becomes p_{2}-p_{1}, in which case the dot product is:

cos(0)*sin(q_{1})*cos(p_{2}-p_{1})*sin(q_{2}) +
sin(0)*sin(q_{1})*sin(p_{2}-p_{1})*sin(q_{2}) +
cos(q_{1})*cos(q_{2}). =

sin(q_{1})*cos(p_{2}-p_{1})*sin(q_{2}) +
cos(q_{1})*cos(q_{2}).

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(q_{1})*cos(p_{2}-p_{1})*sin(q_{2}) +
cos(q_{1})*cos(q_{2}) / (1 * 1). = sin(q_{1})*cos(p_{2}-p_{1})*sin(q_{2}) +
cos(q_{1})*cos(q_{2})].

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(q_{1})*cos(p_{2}-p_{1})*sin(q_{2}) +
cos(q_{1})*cos(q_{2}))] * r.

Michael Shackleford, A.S.A.