From: jeremy@deneb.Princeton.EDU (Jeremy Goodman)
Subject: Re: solving a great circle problem
Originator: news@nimaster
Organization: Princeton University
Distribution: na
Date: Wed, 30 Dec 1992 21:50:34 GMT

(This is rather long and technical.)
Eric Varsanyi writes

>Using the trig in the 1943 book I wrote a great circle distance and
>bearing subroutine, but it seems like there should be an easier
>(less floating point multiply's) way to do this.

>It looks like the fplan program uses a non-great circle routing, I
>would appreciate a pointer to a description of that algorithm if 
>anyone has it as well.

...and he quotes the fplan algorithm as follows:

>
>int
>psolve(m)
>struct measure *m;
>{
>        double tmp,dist,hdg;
>
>   tmp = atan2 (2.0 * asin (sin ((m->to_lon - m->from_lon) / 2.0)),
>                log (tan (PI / 4.0 + m->to_lat / 2.0)) -
>                log (tan (PI / 4.0 + m->from_lat / 2.0)) );
>   if (m->from_lat == m->to_lat)
>      dist = 60.0 *
>         fabs (RAD2DEG (2.0 * asin (sin ((m->to_lon - m->from_lon) / 2.0)))) *
>            cos (m->from_lat);
>   else
>      dist = 60.0 * RAD2DEG (m->to_lat - m->from_lat) / cos (fabs (tmp));
>
>   if (asin (sin (m->to_lon-m->from_lon)) >= 0.0)
>      hdg = RAD2DEG (fabs(tmp));
>   else
>      hdg = 360.0 - RAD2DEG (fabs(tmp));
>
>        printf("(Dist=%10.2lf Hdg=%5.1lf)\n",dist,hdg);
>}
>

The course that fplan is computing is not a great circle, but one
that maintains a constant true heading.  Such a course is called a
"rhumbline" or "loxodrome".  The easiest way to understand such courses
is in terms of a Mercator projection, which maps longitude and latitude
into Cartesian coordinates (x,y) according to the rule

x = longitude
y = log[tan(pi/4 + longitude/2)] 	(natural log here).

This is a conformal map, which means that it preserves angles:
two curves on the sphere intersecting at a given angle are mapped
into curves on the xy plane meeting at the same angle.
Since meridians map into lines parallel to the y axis, it follows that
a rhumbline, which meets all meridians at the same angle, maps into
a straight line in the xy plane; this property no doubt explains why the
Mercator projection has been widely used by sailors ever since it was 
invented in the 16th century!  The meaning of the formula for tmp
(the constant heading) should now be pretty clear, apart perhaps from
	2.0 * asin (sin ((m->to_lon - m->from_lon) / 2.0))
which merely gives the properly signed difference in longitudes along 
the shorter (east or west bound) route (a very compact formula, but both
more obscure and more expensive than an if statement!).

Furthermore, the arc length along the rhumbline is dx/dy=constant

dist = sqrt[1 + (dx/dy)**2]*abs(latitude_2 - latitude_1).

Proving these statements is not too hard using elementary differential 
geometry.  I stunned my instructor during my private training by deriving 
these formulae and programming them into my HP calculator because I found 
it a pain to draw course lines across two sides of a sectional.

Although a constant true heading is convenient, a constant MAGNETIC
heading would be even more so.  To the extent that the earth's magnetic
field can be represented by a simple dipole, this could be achieved by
using a system of latitude and longitude based upon the magnetic (rather
than the rotational) poles, and then carrying out a Mercator projection
from this system.

Eric also asks about efficient computation of the distance c and angle A in the
spherical triangle (L1 to L2 is a great circle)
 *
 *               True North Pole
 *                      +
 *                     /C\
 *                    /   \
 *                 b /     \ a
 *                  /       \
 *                 /A       B\
 *                +-----------+
 *               L1     c      L2
 *
 *      Given:
 *              C - difference in longtitude (shortest)
 *              b - co-latitude of L1 (90-lat)
 *              a - co-latitude of L2 (90-lat)
 *
 *      Find:
 *              c - angular distance between L1 and L2, great circle
 *              A - bearing from L1 to L2, great circle
 *
Eric could save one trig call in the computation of the distance by first
computing

	cplus  = 0.5*cos(a+b)
	cminus = 0.5*cos(a-b).

The cosine formula then reads

	cos(c)= cminus+cplus +(cminus-cplus)*cos(C)

One drawback of the cosine formula is that it will give a relatively 
large roundoff error for very short distances, when cos(c) is
close to unity.  This may not be a problem in practice if double
precision is used.

In the computation of the bearing from L1 to L2, Eric uses the cosine
formula a second time.  With the above method for the distance, it would 
be more efficient to use the sine formula,

sin(A)/sin(a) = sin(C)/sin(c),

since sin(b) will not have been computed.  Still further savings are
possible if one computes sin(C) and cos(C) from

t= tan(C/2); 	sin(C)= 2*t/(1+t*t);	cos(C)= (1-t*t)/(1+t*t);

provided that a single call to tan is less expensive than calling both
sin and cos.
In any case, the heading along a great circle is not constant,
so there will be a need to compute intermediate headings.


I learned my spherical trig from
TEXT-BOOK ON SPHERICAL ASTRONOMY by W.M. Smart (Cambridge Univ Press, 1947),
---which is old-fashioned, but lucid.

A reference on map projections I found in the astronomy library here is
MAP PROJECTION METHODS, by Frederick Pearson (Sigma Scientific, 1984)
---but I am not familiar with this book.  There may be something better.

Finally, I for one would be interested in the final optimized programs
that Eric produces.