Chronology |
Current Month |
Current Thread |
Current Date |

[Year List] [Month List (current year)] | [Date Index] [Thread Index] | [Thread Prev] [Thread Next] | [Date Prev] [Date Next] |

*From*: John Denker <jsd@AV8N.COM>*Date*: Fri, 20 Aug 2004 17:38:59 -0400

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.

- Prev by Date:
**Re: spherical geometry** - Next by Date:
**Re: Expert calls for optional maths** - Previous by thread:
**Free Physics Software** - Next by thread:
**Re: Expert calls for optional maths** - Index(es):