From: ewv@craycos.com (Eric Varsanyi)
Newsgroups: rec.aviation.misc,rec.aviation.ifr
Subject: solving a great circle problem
Date: 26 Dec 92 20:46:06 GMT
Organization: Cray Computer Corporation

I'm working on some GPS moving map software (that runs on an HP95) and
am looking for references to any literature describing spherical
trigonometry. I have a book from 1943 (Navigational Trigonometry,
by Rider and Hutchinson) that lays all the trig groundwork for solving
spherical triangles, but I suspect the state of the art has advanced
since then. A search of various engineering libraries (under navigation+
trigonometry and a few other permutations) keeps popping up this book
and one other (real old) one I can't get my hands on.

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.

I'll post my final (optimized) algorithm to the net if there is interest...

-Eric Varsanyi
 ewv@vars.com

Here is my (un-optimized) great circle solution:
/*
 * Solve the following spherical triangle:
 *
 *               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
 *
 *      West longtitude is negative, South latitude is negative
 */

int
solve(m)
struct measure *m;
{
        double  dlo;            /* Difference in longtitude */
        double  b,a;            /* Co-latitude of lat1 and lat2 */
        double  c;              /* Angular GC ditance between L1 and L2 */
        double  sinb;
        double  brng;
        double  dist,hdg;
        int     wflag = 0;      /* Going west */

        /*
         * Compute DLo, shortest way around
         */

        dlo = fabs(m->from_lon-m->to_lon);
        if (m->to_lon < m->from_lon) wflag = 1;         /* Going west */
        if (dlo > DEG2RAD(180.0)) {
                /* Quicker to go the other way */
                dlo = DEG2RAD(360.0) - dlo;
                wflag = 1-wflag;
        }

        /*
         * Compute co-latitudes
         */
        b = DEG2RAD(90.0)-m->from_lat;
        a = DEG2RAD(90.0)-m->to_lat;

        /*
         * Compute angular distance
         */
        sinb = sin(b);          /* Used twice */
        c = acos(cos(b-a)-(sin(a)*sinb*(1-cos(dlo))));

        /*
         * Convert to NM
         */
        dist = c*(60*R2D);      /* 1 minute == 1 NM */

        /*
         * Compute bearing from L1 to L2
         */
        brng = acos( ((cos(a)-cos(c-b))/(sinb*sin(c)))+1 );

        /*
         * Normalize if we're going west
         */
        if (wflag) brng = DEG2RAD(360.0) - brng;

        /*
         * Compute to degrees heading
         */
        hdg = RAD2DEG(brng);

        printf("%s->%s Dist=%10.2lf Hdg=%5.1lf ",m->from,m->to,dist,hdg);
}

Here's the 'fplan' algorithm, I understand the asin(sin) stuff (normalizing
the inputs) but the 'tmp = atan2(xxx)' is a mystery.

/*
 * Solve using (projected rectangular co-ords?)
 */
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);
}