Here is my solution to the great-circle navigation problem.
This problem is a veritable poster-child for illustrating
the power and elegance of the Clifford Algebra approach
... which could also be called the quaternion approach,
since quaternions are isomorphic to unit bivectors in
Clifford Algebra. So if you can do bivectors, you get
quaternions for free.
Statement of the problem: We are given two points on the
sphere, a and b. The task is to find the heading that we
should use when departing a enroute to b.
WLoG we scale things so that our sphere is the unit sphere.
Consequently a and b are unit vectors.
Assumption #1: We assume a is not a multiple of b, because
otherwise the problem is trivial ... *any* direction will
equally well get you from a to b. To say the same thing
more formally, we assume a/\b is nonzero.
Assumption #2: We assume a is not the north or south pole.
Otherwise expressing the answer as a heading is not useful.
(Some useful alternatives will be given below.)
By way of preparation, it will be convenient to expand b
as follows:
b := c a + s w
where w is some vector that lies in the a/\b plane and
is perpendicular to a, i.e. w.a = 0. Here c and s
are scalars. (They are the cosine and sine of some
angle, but that's actually more than we need to know.)
You can easily verify that a/\b = s a/\w. That makes
sense; the a/\b plane is the same as the a/\w plane,
except for a scalar scale-factor.
Thence you can easily verify that a(a/\b) = s w, or
equivalently a.(a/\b) = s w. That gives us a
convenient way of constructing w:
w = a.(a/\b) / ||a/\b|| [eq 1]
... where the denominator is well-behaved by virtue
of assumption #1 (above).
As usual the norm is defined as
||a/\b|| := sqrt((a/\b).(b/\a)) [eq 2]
So much for preparation. Let's actually attack the
problem now.
Let x(theta) be a point that roves along the great
circle from a to b, starting with x=a when theta=0.
Then the tangent vector will be d(x)/d(theta).
My physical intuition tells me that the tangent vector,
tangent to the point a, will be just w (or perhaps a
multiple thereof, depending on how you scale theta).
(If that's not obvious to you, see the appendix, below.)
So here's the recipe: If the two vectors a and b are
already expressed in rectangular coordinates, that's
great. Otherwise, e.g. if they are in polar coordinates
(latitude and longitude), convert them to rectangular
coordinates ... which will be the last (or almost the
last) use of trig functions in the whole problem!
Find the vector w using equation [1] as given above.
Next we need to know what w looks like to the guy on
the ground, i.e. in the tangent plane, tangent to the
sphere at point a. To get our bearings, let's determine
the vector (call it v) that points from a to the north
pole. We can find this immediately by using equation
[1] again ... just imagine traveling from a to the north
pole.
Now we project out the component of w parallel to v (call it
w1). The remainder of w (call it w2) will be perpendicular
to v. All of these vectors (w, v, w1, and w2) lie in the
tangent plane. We use the usual formula for the projection
operation:
w1 = v (v.w) / v.v
and
w2 = w - w1
At this point we know enough to draw the vector that points
from a to b. Lay a piece of graph paper on the floor
(aligned NS/EW of course), and draw a vector from the origin
to the point (w2,w1).
Note that starting from the rectangular components of a
and b, we have gotten to this point without using any trig
functions or any other transcendental functions. We have
just performed algebraic computations on the components.
It's just add, subtract, multiply, and divide. We can
even leave out the square roots, since we only care about
the direction of w, not its magnitude. This is arguably
important if you need to do this calculation a zillion
times per second, and you don't want to spend all your CPU
time evaluating transcendental functions.
In general there are two great-circle routes from a to b:
a short one and a long one. The formalism given here is
not guaranteed to give you the short one. There is an
ambiguity associated (conceptually, at least) with the
choice of positive versus negative square root in
equation [2]. Or you could just put an explicit +-
sign in the definition of w in equation [1]. (I suspect
it's possible to set things up so that the formalism
always coughs up the short route, but I haven't pursued
this idea.)
Finally: If you insist on expressing the answer as a
heading angle, calculate atan2(w1,w2).
==========
For slight extra credit, we now consider the special case
where a is at the north pole. The desired departure heading
is south (duh!) but that probably isn't what you wanted to
know. To determine which meridian you should follow, you
can skip the entire calculation; the answer is just the
longitude of b, by inspection.
However, you can appreciate the power and beauty of the
bivector formulation by observing that you *can*
calculate the vector w just fine, even when a is at the
north pole. We cannot calculate the projections w1 and
w2, because v doesn't exist ... but w exists as a
perfectly well-behaved vector in real space.
This is important if you're building something like an
autopilot or video game. The flight director can display
the vector w as a pointer that points that-a-way ... and
it will work just fine at the poles and everywhere else.
To say the same thing in fancier terms:
-- expressions involving polar coordinates are singular
at the poles and numerically unstable near the poles, but
-- the actual physics and geometry of the sphere have
no singularities anywhere, and
-- the bivector representation has no singularities
anywhere.
============================================
Appendix.
You can easily prove that the tangent to the great
circle at point a is just w, a vector in the a/\b
plane perpendicular to a.
This should be obvious from the geometry of the situation,
but it is also possible to prove it by directly calculating
the tangent vector. It's an amusing exercise in Clifford
Algebra. We express the roving point x(theta) by applying
the rotation operator R(theta) to the vector a, and then
differentiating w.r.t theta. Of course we are talking about
the particular rotation operator R() that causes rotations
in the a/\b plane.
Specifically, for an infinitesimal angle theta, we can
use the standard formula http://www.av8n.com/physics/rotations.htm#eq-rotorize-gamma
for producing a rotation, which gives us:
R(theta) a = [1 + (theta/2)w/\a] a [1 + (theta/2)a/\w]
and if you just turn the crank, expanding the RHS and
simplifying using the axiomatic properties of vectors
a and w, you get
... = a + theta w + irrelevant terms involving theta^2
This can be immediately differentiated w.r.t theta,
whereupon we find that the tangent vector is w, QED.