The history of geodesic calculators at APL/Uw (and beyond)
is a strange odyssey. A geodesic calculator, which outputs 
a geodesic track on, in this case, an oblate spheroid model
of the Earth, is a fundamental building block of long range
ocean acoustic modeling packages. The beginnings of the
effort that lead to the current package are obscured in the
mists of time, but probably involved Peter Worcester at Scripps  
during the inception of long range ocean acoustic tomography.
Fossilized relics of that time persist in comment code fragments
authored by Bruce Cornuelle (accusing Peter of something.)
Recently further evidence has come to light to this author
(on a climb of Mt. Shuksan in the N. Cascades) that suggests
that Cornuelle wrote a Fortran routine that performed geodesic 
calculations. Knowing Bruce, that would be about right.

This Fortran code was passed at some time to Rich Pawlowicz 
(the other person on the Shuksan climb with me....) who was
a graduate student at WHOI at the time. Rich ported the 
Fortran to Matlab, embellishing a little, and so on.

The Worcester-Cornuelle code was known in the early 90's 
to APL/UW when we undertook to build our own long range 
propagation modeling package. This was the original MAP
program. [S. Leach, OCEANS something]. George Dworski, a 
mathematician at APL/UW reviewed the Worcester-Cornuelle code
and reimplemented it in Fortran, using an alternate numerical
approach. (More on this later.) Dworski also included the 
ability for the calculations to include a data base of 
ocean sound speed profiles, and thereby compute refracted
geodesic paths. 

Dworski's code was cleaned up and polished by graduate student
Craig Friesen, and the resulting executables formed the geodesic
calculator for the first APL/UW MAP program.

It was then decided at APL/UW to reimplement the MAP program
in Matlab [C. Eggen, OCEANS something]. The Matlab version
used the Pawlowicz routine as-is. It still remains in the state
that Rich left, oh so many years ago.

In about 2000, undergraduate student Mark Peeples and I revisited
the MAP program (in its Matlab reincarnation) and determined that
we should port parts of it out of Matlab language. Mark undertook
the laborious and thankless effort of porting the Pawlowicz
code into C.  This was a difficult effort, but Mark produced a 
routine that entered first level unit testing. We tested it against
the Matlab routine for a number of test cases and discovered to
our surprise that the Matlab routine simply bombed for a number of
simple (or so we thought) cases. 

And then I remembered that conversation on the side of Mt. Shuksan
where Rich recalled that his routine occasionally had problems of
various sorts that he could never quite debug....

At this point I revisited the two approaches (the Worcester-Cornuelle-
Pawlowicz or WCP approach and the Dworski-Friesen of DF approach) and 
made some interesting observations, to wit:

1) The WCP code implements a hybrid of algorithms found in the tome
Geodesy by Bomgren. This book, originally published in something like
1910, is the first word in geodesy. It devotes at least one huge chapter on 
a variety of approximations to solving the geodesic problem. There
is a brief mention that the equations might be solved by digital computers
that bigger institutions may have at their disposal, and outlines a
way to do this. The bulk of the algorithms, however, are designed
for workers in the field who do not have high-powered numerical
solvers, and therefore have to rely on mathematically simpler approximations
valid to various orders. It is these approximations that the 
WCP code implements.

2) The DF approach tackles the geodesic problem directly. The 1-point
problem (given a point, a bearing and a range, find the other point)
is a shooting problem, and the 2-point problem (given two points, 
find the range of the shortest path between, and the forward and reverse
bearings at either end) is the brachistochrone problem in calculus
of variations. However, on an oblate spheroid, after some elegant
mathematics, the variational problem can be reduced to two coupled
ordinary differential equations. This was Dworski's insight: that
both problems would be "solved" using a shooting method solution
to ODEs.

The difference between the two approaches is this: the WCP code
implements an approximate solution exactly, whereas the DF code
implements the exact solution with numerical approximations that
have well-understood convergence properties.

There was a drawback to the DF code, in that it did not handle
meridional problems at all. Fortran is also not entirely portable,
so I took the work that Mark and I started and developed a full-fledged
object-oriented C++ class that solves the 1 and 2 point problems
using the same techniques as the DF code. In short, the package
provides a geodetic calculator implemented as a C++ class which
compiles into a library. The class is based on a Runge-Kutta integrator
which solves the shooting problem directly. The brachistochrone
problem uses a Nelder-Mead simplex minimizer that contains the
Runge-Kutta integrator and basically shoots over a solution space
of range and bearing until the "shot" hits the second point.

Our experience with problems and testing inspired a suite of test
codes. There are 2 test generators in the package, gen_shoot and 
gen_brach, which generate random problems for the 1-point and 
2-point problems, respectively.  There are 2 test harnesses, 
test_shoot and test_brach, which take in the test problems and 
apply the geodesic calculator to produce output for the 
1-point and 2-point problems, respectively. And there are 2
comparators, cmp_shoot and cmp_brach, which each take in a pair
of outputs (from whatever source) and output the errors between
calculations.

(I recently added a kind of stand-alone tester, test_shoot_recip/,
which shoots from P1 to P2, and then back to P1 (to, say, P3) and
then determines how close P3 was to P1. The Matlab dist() routine
was known to have some problems doing this kind of reciprocal problem,
Pawlowicz, personal communication I think.)

The geoids implemented by the class are: spherical, WGS84, 
Clarke66, IAU73. 

The class also has fully implemented a solution involving
oblique spherical triangles for a spherical earth
for both the 1-point and 2-point problem.

Both types of solutions (integrated geodesic equations and spherical
triangles) are supposed to handle all cases, including singular, 
meridional, polar and regular non-meridional. There isn't a case
that isn't supposed to be handled somewhere. In general, the code
to handle all the special or weird cases (usually involving a branch
crossing) increased the code size by about 300%. Extra error handling
and reporting and associated class administrivia routines add
another 50% to the class code base.

There is also a directory of basic test harnesses for each case.

Testing: Testing is a primary concern, and the philosophy here is
to provide both a standard set of results and the tools to test
this or other codes against each other or against the standard
results. The first test that can be run is the calculator in
integrator mode versus the calculator in spherical triangle mode
for a spherical earth. Here, the solution from the spherical
triangle code is taken to be true, and errors in comparison are
attributed to the numerical Runge-Kutta and Nelder-Mead routines.
The integrator with its standard "default" tolerances 
holds up very well against the spherical triangle
solution "standard" for all test cases for both the 1-point and 
2-point problems: the summary is in the directory "test_cases". 
For the 1-point problem, the solutions for final latitude,
longitude and reverse bearing (at final point) have errors of
about 1e-4 degrees or less (often orders of magnitude less.)
Similarly for the 2-point case, the forward and reverse bearings
have errors of about 1e-4 degrees or less, and the range errors
are about 1 - 5 metres. 

There is no standard result for geodesics computed on an oblate
spheroid, so the approach taken here is to use the solution provided
by the class, after it has been shown to perform well on a spherical
earth, and call that the standard. Errors between this "standard" 
and the true values should be approximately the same order of magnitude
as between the integrator and the spherical trig results for the
spherical earth.

I have compared the WGS84 solution provided by the class against
that provided by Dworski's code, for the 2-point problem. Since
both codes use the same mathematical algorithms, this test primarily 
compares numerical implementations of these algorithms. (Note also
that Dworski's code does not handle meridional cases.) Since the
2-point problem is solved using minimizers and shooters, this
comparison more or less compares both the 1-point (shooter) code
and the 2-point code. The errors in range were O(10m) or less, and
the errors in forward bearing were O(0.001 deg) or less. (Dworski's
code does not report final reverse bearing.)

This comparison suggests that this class is faithful to Dworski's
solution for non-meridional 1-point and 2-point problems on an
oblate spheroid. 

The comparison between the integrator and the WCP Matlab routine 
"dist.m", which the current MAP program is based on, is interesting:

1-point: The "dist" routine does not do shooting problems, so 
no comparison is possible.

2-point: Two comparisons were made: for a spherical earth
and the WGS84 spheroid:

2-point::sphere: dist had problems with some meridional cases and/or
cases involving one or the other pole. This behaviour was known
from the testing Mark and I did. For non-meridional cases, 
the errors were comparable to those between the integrator and 
the spherical trig solutions.

2-point::wgs84: Forward and reverse bearing estimates were usually
good to within 0.5 degrees or better, sometimes much better. There
was a considerable range in range errors, ranging from 1 metre to
as much as 200km! There were some severe anomalies associated with 
problems involving points on poles, as noted above. I did not try to
quantify which classes of 2-point problems had severe range errors
and which did not.

The "gold standard" for oblate spheroids: This testing does not say
whether the DF solution or the WCP solution is more accurate for
oblate spheroids. When the "oblate spheroid" is reduced to the simple
case of a perfect sphere, both solutions work well. There is an inherent
elegance to the DF solution in that it solves the exact problem with
numerical codes that can be adjusted to have ever tighter error tolerances.
In principle, the code could be rerun with ever tighter tolerance
parameters and thereby point to a "convergent" solution which must
surely be the right solution.

The WCP code, based on approximations, cannot do this.

I have not pursued such a convergence study. There are methods in the class that
expose the integrator and minimizer tolerance parameters, so one could 
pursue such a task. (Give it to an undergrad, Mark Peeples, where are you?) 

It should also be noted that the geodesy calculator handled all the meridional
polar and non-polar cases, whereas the DF code did not, and the WCP
occasionally had some wild errors. 


The results from testing the dist.m routine are also stored in the
test_cases directory.