[Proj] Accurate algorithm for geodesic calculations

[Charles Karney]

I've implemented an algorithm for geodesic calculations for an ellipsoid
with the following features:

 * It uses expansions accurate to 8th order in the ellipsoid flattening.
   This gives 12 nm accuracy using doubles and 6 pm accuracy using long
   doubles.
 * The direct calculation uses a reverted series and is therefore
   non-iterative.  A mechanism is provided to do a series of direct
   calculations on a single geodesic which is faster by a factor of 2.5.
 * The inverse calculation uses Newton's method to determine the
   azimuth.  This typically converges in 3 iterations.  I haven't found
   any cases where it fails to converge.
 * I provide a large test set geodesics for the WGS84 ellipsoid which is
   accurate to 0.1 pm and 10^-18 degree.  Thus was generated by Maxima
   using 20th order expansions (approximately accurate to 1 part in
   10^50).

Further information is available at

    http://charles.karney.info/geographic/geodesic.html

The code itself (C++ classes and a program for doing geodesic
calculations) is available at

    http://charles.karney.info/geographic

This compiles with g++ and Visual Studio 2005.

Still to do:

 * Reduce the order of the expansions to match precision available with
   doubles.
 * Improve the starting guesses for Newton's method to reduce to number
   of iterations in some cases.
 * More documentation on the method.
 * Provide Maxima code to generate the expansions.

-- 
Charles Karney <ckarney at sarnoff.com>
Sarnoff Corporation, Princeton, NJ 08543-5300

URL: http://charles.karney.info
Tel: +1 609 734 2312
Fax: +1 609 734 2662