// -*- text -*-
/**
 * \file Geographic.dox
 * \brief Documentation for GeographicLib
 *
 * Written by Charles Karney <charles@karney.com> and licensed under the
 * MIT/X11 License.  For more information, see
 * http://geographiclib.sourceforge.net/
 **********************************************************************/

/**
\mainpage Geographic library
\author Charles F. F. Karney (charles@karney.com)
\version 1.30
\date 2013-02-27

\section abstract Abstract

%GeographicLib is a small set of <a href="annotated.html">C++
classes</a> for performing conversions between geographic, UTM, UPS,
MGRS, geocentric, and local cartesian coordinates, for gravity (e.g.,
EGM2008), geoid height and geomagnetic field (e.g., WMM2010)
calculations, and for solving geodesic problems.  The emphasis is on
returning accurate results with errors close to round-off (about
5--15 nanometers).  New accurate algorithms for \ref geodesic and
\ref transversemercator have been developed for this library.  The
functionality of the library can be accessed from user code, from the
\ref utilities provided, or via the \ref other.  For a sample of the
geodesic capabilities in Javascript, check out the
<a href="../scripts/geod-calc.html">online geodesic calculator</a> and
the script for displaying
<a href="../scripts/geod-google.html">geodesics in Google Maps</a>

This library is <i>not</i> a general purpose projection library; use
<a href="http://trac.osgeo.org/proj/">proj.4</a> for that.  On the other
hand, it does provide the core functionality offered by
<a href="http://earth-info.nima.mil/GandG/geotrans/">GEOTRANS</a>.

\section download Download

The main project page is at
- <a href="http://sourceforge.net/projects/geographiclib">
  http://sourceforge.net/projects/geographiclib </a>
.
The code is available for download at
- <a href="http://sf.net/projects/geographiclib/files/distrib/GeographicLib-1.30.tar.gz">
  GeographicLib-1.30.tar.gz</a>
- <a href="http://sf.net/projects/geographiclib/files/distrib/GeographicLib-1.30.zip">
  GeographicLib-1.30.zip</a>
.
as either a compressed tar file (tar.gz) or a zip file.  (The two
archives have identical contents, except that the zip file has DOS
line endings.)  Alternatively you can get the latest release using git
\verbatim
  git clone git://git.code.sf.net/p/geographiclib/code geographiclib
  cd geographiclib; git checkout -q r1.30
\endverbatim
There are also binary installers for Windows available at
- <a href="http://sf.net/projects/geographiclib/files/distrib/GeographicLib-1.30-win32.exe">
  GeographicLib-1.30-win32.exe</a>
- <a href="http://sf.net/projects/geographiclib/files/distrib/GeographicLib-1.30-win64.exe">
  GeographicLib-1.30-win64.exe</a>
.
It is licensed under the
<a href="http://www.opensource.org/licenses/MIT">MIT/X11 License</a>;
see <a href="LICENSE.txt">LICENSE.txt</a> for the terms.
For more information, see http://geographiclib.sourceforge.net/.

\section contents Contents
 - \ref intro
 - \ref install
 - \ref start
 - \ref utilities
 - \ref organization
 - \ref other
 - \ref geoid
 - \ref gravity
 - \ref magnetic
 - \ref geodesic
 - \ref triaxial
 - \ref transversemercator
 - \ref geocentric
 - \ref old

<center>
Forward to \ref intro.
</center>

**********************************************************************/
/**
\page intro Introduction

<center>
Forward to \ref install.  Up to \ref contents.
</center>

%GeographicLib offers a C++ interfaces to a small (but important!) set
of geographic transformations.  It grew out of a desire to improve on
the <a href="http://earth-info.nima.mil/GandG/geotrans/">geotrans</a>
package for transforming between geographic and MGRS coordinates.  At
present, %GeographicLib provides UTM, UPS, MGRS, geocentric, and local
cartesian projections, gravity and geomagnetic models, and classes for
geodesic calculations.

The goals of %GeographicLib are:
 - Accuracy.  In most applications the accuracy is close to round-off,
   about 5 nm (5 nanometers).  Even though in many geographic
   applications 1 cm is considered "accurate enough", there is little
   penalty in providing much better accuracy.  In situations where a
   faster approximate algorithm is necessary, %GeographicLib offers an
   accurate benchmark to guide the development.
 - Completeness.  For each of the projections included, an attempt is
   made to provide a complete solution.  For example,
   GeographicLib::Geodesic::Inverse works for anti-podal points.
   Similarly, GeographicLib::Geocentric.Reverse will return accurate
   geodetic coordinates even for points close to the center of the
   earth.
 - C++ interface.  For the projection methods, this allows encapsulation
   of the ellipsoid parameters.
 - Emphasis on projections necessary for analyzing military data.
 - Uniform treatment of UTM/UPS.  The GeographicLib::UTMUPS class treats
   UPS as zone 0.  This simplifies conversions between UTM and UPS
   coordinates, etc.
 - Well defined and stable conventions for the conversion between
   UTM/UPS to MGRS coordinates.
 - Detailed internal documentation on the algorithms.  For the most part
   %GeographicLib uses published algorithms and references are given.  If
   changes have been made (usually to improve the numerical accuracy),
   these are described in the code.

Various \ref utilities are provided with the library.  These illustrate
the use of the library and are useful in their own right.  This library
and the utilities have been tested with g++ 4.4 under Linux, with g++
4.2 under Mac OS X, and with MS Visual Studio 9 (2008), 10 (2010), and
11 (2012) compiled for 32 bit and 64 bit.  Note that one of the tests \b
fails with Visual Studio 11.  This is a bug in the I/O library and is
unlikely to affect most users (bug reported 2013-01-10).

Matlab, Javascript, and Python interfaces are provided to portions of
%GeographicLib; see \ref other.

The section \ref geodesic documents the method of solving the geodesic
problem.

The section \ref transversemercator documents various properties of this
projection.

The bulk of the testing has used geographically relevant values of the
flattening.  Thus, you can expect close to full accuracy for -0.01 &le;
\e f &le; 0.01 (but note that GeographicLib::TransverseMercatorExact is
restricted to \e f > 0).  However, reasonably accurate results can be
expected if -0.1 &le; \e f &le; 0.1.  Outside this range, you should
attempt to verify the accuracy of the routines independently.  Two types
of problems may occur with larger values of \e f:
 - Some classes, specifically GeographicLib::Geodesic,
   GeographicLib::GeodesicLine, and GeographicLib::TransverseMercator,
   use series expansions using \e f as a small parameter.  The accuracy
   of these routines will degrade as \e f becomes large.
 - Even when exact formulas are used, many of the classes need to invert
   the exact formulas (e.g., to invert a projection), typically, using
   Newton's method.  This usually provides an essentially exact
   inversion.  However, the choice of starting guess and the exit
   conditions have been tuned to cover small values of \e f and the
   inversion may be incorrect if \e f is large.

Undoubtedly, bugs lurk in this code and in the documentation.  Please
report any you find to charles@karney.com.

<center>
Forward to \ref install.  Up to \ref contents.
</center>
**********************************************************************/
/**
\page install Installing %GeographicLib

<center>
Back to \ref intro.  Forward to \ref start.  Up to \ref contents.
</center>

%GeographicLib has been developed under Linux with the g++ compiler (versions
4.0 and later) and under Windows with Visual Studio 2005, 2008, and
2010.  Earlier versions were tested also under Darwin and Solaris.  It
should compile on a wide range of other systems.  First
download either
<a href="http://sourceforge.net/projects/geographiclib/files/distrib/GeographicLib-1.30.tar.gz">
GeographicLib-1.30.tar.gz</a> or
<a href="http://sourceforge.net/projects/geographiclib/files/distrib/GeographicLib-1.30.zip">
GeographicLib-1.30.zip</a> (or
<a href="http://sourceforge.net/projects/geographiclib/files/distrib/GeographicLib-1.30-win32.exe">
GeographicLib-1.30-win32.exe</a> or
<a href="http://sourceforge.net/projects/geographiclib/files/distrib/GeographicLib-1.30-win64.exe">
GeographicLib-1.30-win64.exe</a> for binary installation under Windows).
Then pick one of the first five options below:
- \ref cmake.  This is the preferred installation method as it will work
  on the widest range of platforms.  However it requires that you have
  <a href="http://www.cmake.org">cmake</a> installed.
- \ref autoconf.  This method works for most Unix-like systems,
  including Linux and Mac OS X.
- \ref gnu.  This is a simple installation method that works with g++
  and GNU make on Linux and many Unix platforms.
- \ref windows.  This is a simple installation method that works with
  Visual Studio 2005, 2008, and 2010 under Windows.
- \ref windowsbin.  Use this installation method if you only need to use
  the \ref utilities supplied with %GeographicLib.  (This method also
  installs the header files and the static library for use by Visual
  Studio 2010.)
- \ref qt.  How to compile %GeographicLib so that it can be used by Qt
  programs.
- \ref maintainer.  This describes addition tasks of interest only to
  the maintainers of this code.
.
This section documents only how to install the <i>software</i>.  If you
wish to use %GeographicLib to evaluate geoid heights or the earth's
gravitational or magnetic fields, then you must also install the
relevant data files.  See \ref geoidinst, \ref gravityinst, and \ref
magneticinst for instructions.

The first two installation methods use two important techniques which
make software maintanence simpler
- <b>Out-of-source builds:</b>  This means that you create a separate
  directory for compiling the code.  In the description here the
  directories are called BUILD and are located in the top-level of the
  source tree.  You might want to use a suffix to denote the type of
  build, e.g., BUILD-vc9 for Visual Studio 9, or BUILD-shared for a
  build which creates a shared library.  The advantages of out-of-source
  builds are:
  - You don't mess up the source tree, so it's easy to "clean up".
    Indeed the source tree might be on a read-only file system.
  - Builds for multiple platforms or compilers don't interfere with each
    other.
- <b>The library is installed:</b> After compilation, there is a
  separate <i>install</i> step which copies the headers, libraries,
  tools, and documentation to a "central" location.  You may at this
  point delete the source and build directories.  If you have
  administrative privileges, you can install %GeographicLib for the use
  of all users (e.g., in /usr/local).  Otherwise, you can install it for
  your personal use (e.g., in $HOME/packages).

\section cmake Installation with cmake

This is the recommended method of installation; however it requires that
<a href="http://www.cmake.org">cmake</a> be installed on your system.
This permits %GeographicLib to be built either as a shared or a static
library on a wide variety of systems.  cmake can also determine the
capabilities of your system and adjust the compilation of the
libraries and examples appropriately.

cmake is available for most computer platforms.  On Linux systems cmake
will typically one of the standard packages and can be installed by a
command like
  \verbatim
  yum install cmake \endverbatim
(executed as root).  On other systems, download a binary
installer from http://www.cmake.org click on download, and save and run
the appropriate installer.  Run the cmake command with no arguments to
get help.  Other useful tools are ccmake and cmake-gui which offer
curses and graphical interfaces to cmake.  Building under cmake depends
on whether it is targeting an IDE (interactive development environment)
or generating Unix-style makefiles.  The instructions below have been
tested with makefiles and g++ on Linux and with the Visual Studio IDE on
Windows.

Here are the steps to compile and install %GeographicLib:
- Unpack the source, running one of \verbatim
  tar xfpz GeographicLib-1.30.tar.gz
  unzip -q GeographicLib-1.30.zip \endverbatim
  then enter the directory created with one of \verbatim
  cd GeographicLib-1.30 \endverbatim
- Create a separate build directory and enter it, for example, \verbatim
  mkdir BUILD
  cd BUILD \endverbatim
- Run cmake, pointing it to the source directory (..).  On Linux, Unix,
  and MacOSX systems, the command is \verbatim
  cmake .. \endverbatim
  For Windows, the command is typically one of \verbatim
  cmake -G "Visual Studio 10" -D CMAKE_INSTALL_PREFIX=C:/pkg-vc10/GeographicLib-1.30 ..
  cmake -G "Visual Studio 9 2008" -D CMAKE_INSTALL_PREFIX=C:/pkg-vc9/GeographicLib-1.30 ..
\endverbatim
  Instead of setting CMAKE_INSTALL_PREFIX, you can (on Windows) set
  PACKAGE_PATH to point to its parent directory, e.g., \verbatim
  cmake -G "Visual Studio 10" -D PACKAGE_PATH=C:/pkg-vc10 ..
\endverbatim
  The definitions of CMAKE_INSTALL_PREFIX are optional (see below).  The
  settings given above are recommended to ensure that packages that use
  %GeographicLib use the version compiled with the right compiler.
  If you need to rerun cmake, use \verbatim
  cmake . \endverbatim
  possibly including some options via <code>-D</code> (see the next step).
- cmake allows you to configure how %GeographicLib is built and installed by
  supplying options, for example \verbatim
  cmake -D CMAKE_INSTALL_PREFIX=/tmp/geographic . \endverbatim
  The options you might need to change are
  - <code>COMMON_INSTALL_PATH</code> governs the installation
    convention.  If it is on ON (the Linux default), the installation
    is to a common directory, e.g., /usr/local.  If it is OFF (the
    Windows default), the installation directory contains the package
    name, e.g., C:/pkg/GeographicLib-1.30.  The installation
    directories for the documentation, cmake configuration, python and
    matlab interfaces all depend on the variable with deeper paths
    relative to CMAKE_INSTALL_PREFIX being used when it's ON:
    - documentation: OFF: doc/html; ON: share/doc/GeographicLib/html;
    - cmake configuration: OFF cmake; ON: share/cmake/GeographicLib-1.30;
    - python interface: OFF: python; ON: lib/python/site-packages;
    - matlab interface: OFF: matlab; ON: libexec/GeographicLib/matlab;
    .
  - <code>CMAKE_INSTALL_PREFIX</code> (default: /usr/local on
    non-Windows systems, C:/Program Files/GeographicLib on Windows
    systems) specifies where the library will be installed.  For windows
    systems, it is recommended to use a prefix which includes the
    compiler version, as shown above (and also, possibly, whether this
    is a 64-bit build, e.g., <code>cmake -G "Visual Studio 10 Win64" -D
    CMAKE_INSTALL_PREFIX=C:/pkg-vc10-x64/GeographicLib-1.30 ..</code>).
    If you just want to try the library to see if it suits your needs,
    pick, for example, <code>CMAKE_INSTALL_PREFIX</code>=/tmp/geographic.
  - <code>PACKAGE_PATH</code> and <code>INSTALL_PATH</code> provide an
    alternative mechanism for setting the find_package search path and
    the installation directory.  (find_package is not used by
    %GeographicLib since it doesn't depend on other packages.  However
    PACKAGE_PATH is used here for uniformity with other packages which
    adopt the same conventions.)
    - If PACKAGE_PATH is defined, it is prepended to CMAKE_PREFIX_PATH.
    - If INSTALL_PATH is not specified but PACKAGE_PATH is, then
      INSTALL_PATH is set to
      ${PACKAGE_PATH}, if COMMON_INSTALL_PATH is ON, and to
      ${PACKAGE_PATH}/GeographicLib-1.30, otherwise.
    - If INSTALL_PATH is now defined, then set CMAKE_INSTALL_PREFIX to
      INSTALL_PATH.
    .
    Typically, only PACKAGE_PATH needs to be specified, e.g.,
    - Linux: cmake -D PACKAGE_PATH=/opt .. => CMAKE_INSTALL_PREFIX=/opt
    - Windows: cmake -D PACKAGE_PATH=C:/pkg .. =>
      CMAKE_INSTALL_PREFIX=C:/pkg/GeographicLib-1.30
    .
  - <code>GEOGRAPHICLIB_DATA</code> (default:
    /usr/local/share/GeographicLib for non-Windows systems,
    C:/Documents and Settings/All Users/Application
    Data/GeographicLib for Windows systems) specifies the default
    location for the various datasets for use by GeographicLib::Geoid,
    GeographicLib::GravityModel, and GeographicLib::MagneticModel.
    See \ref geoidinst, \ref gravityinst, and \ref magneticinst for more
    information.
  - <code>GEOGRAPHIC_SHARED_LIB</code> (default: ON for non-Windows
    systems, OFF for Windows systems) determines whether a shared or
    static library is built.  The default is OFF for Windows because
    there's no standard place to install the dlls.
  - <code>CMAKE_BUILD_TYPE</code> (default: Release).  This
    flags only affects non-IDE compile environments (like make + g++).
    The default is actually blank, but this is treated as
    Release.  Choose one of
    \verbatim
  Debug
  Release
  RelWithDebInfo
  MinSizeRel
\endverbatim
    (With IDE compile environments, you get to select the build type in
    the IDE.)
  - <code>MATLAB_COMPILER</code> (default: OFF).  Set this to either
    "mex" (for Matlab) or "mkoctfile" (for Octave) to specify the
    compiler to use for the Matlab/Octave interface.  See \ref matlab
    for more information.
  - <code>GEOGRAPHICLIB_EXAMPLES</code> (default: ON for non-Windows
    systems and OFF for Windows systems).  If set to ON, then the
    short example programs in the <code>examples</code> directory are
    built.
  - <code>GEOGRAPHICLIB_DOCUMENTATION</code> (default: ON for
    non-Windows systems and OFF for Windows systems).  If set to ON,
    then html documentation is created from the source files, provided a
    sufficiently recent version of doxygen can be found.  Otherwise, the
    html documentation will redirect to the appropriate version of the
    online documentation.
  - <code>GEOGRAPHICLIB_PRECISION</code> specifies the precision to be
    used for "real" (i.e., floating point) numbers.  1 means float
    (single precision); 2 (the default) means double; 3 means long
    double; 4 is reserved for quadruple precision.  Nearly all the
    testing has been carried out with doubles and that's the
    recommended configuration.  In particular you should avoid
    "installing" the library with a precision different from double.
- Build and install the software.  In non-IDE environments, run
  \verbatim
  make         # compile the library and the examples
  make test    # run some tests
  make install # as root, if CMAKE_INSTALL_PREFIX is a system directory
\endverbatim
  On IDE environments, run your IDE (e.g., Visual Studio), load
  GeographicLib.sln, pick the build type (e.g., Release), and select
  "Build Solution".  If this succeeds, select "RUN_TESTS" to build;
  finally, select "INSTALL" to install (RUN_TESTS and INSTALL are in
  the CMakePredefinedTargets folder).  Alternatively, you run the
  Visual Studio compiler from the command line with \verbatim
  cmake --build . --config Release --target ALL_BUILD
  cmake --build . --config Release --target RUN_TESTS
  cmake --build . --config Release --target INSTALL
\endverbatim
- The headers, library, and utilities are installed in the
  include/GeographicLib, lib, and bin directories under
  <code>CMAKE_INSTALL_PREFIX</code>.  (dll dynamic libraries are
  installed in bin.)  For documentation, open in a web browser
  <a href="index.html">
  PREFIX/share/doc/GeographicLib/html/index.html</a>, if
  COMMON_INSTALL_PATH is ON, or <a href="index.html">
  PREFIX/doc/index.html</a>, otherwise.

\section autoconf Installation using the autoconfigure tools

The method works on most Unix-like systems including Linux and Mac OS X.
Here are the steps to compile and install %GeographicLib:
- Unpack the source, running \verbatim
  tar xfpz GeographicLib-1.30.tar.gz \endverbatim
  then enter the directory created \verbatim
  cd GeographicLib-1.30 \endverbatim
- Create a separate build directory and enter it, for example, \verbatim
  mkdir BUILD
  cd BUILD \endverbatim
- Configure the software, specifing the path of the source directory,
  with \verbatim
  ../configure \endverbatim
- By default %GeographicLib will be installed under /usr/local.
  You can change this with, for example \verbatim
  ../configure --prefix=/tmp/geographic \endverbatim
- Compile and install the software with \verbatim
  make
  make install \endverbatim
- The headers, library, and utilities are installed in the
  include/GeographicLib, lib, and bin directories under
  <code>prefix</code>.  This installation method does not compile
  the Matlab/Octave interface;  however the source for the interface is
  installed in libexec/GeographicLib/matlab, see \ref matlab of
  instructions on compiling the interface.  For documentation, open
  <a href="index.html">
  share/doc/GeographicLib/html/index.html</a> in a web browser.

\section gnu Installation with GNU compiler and Make

This method requires the standard GNU suite of tools, in particular make
and g++.  This builds a static library and the examples.

Here are the steps to compile and install %GeographicLib:
- Unpack the source, running \verbatim
  tar xfpz GeographicLib-1.30.tar.gz \endverbatim
  then enter the directory created \verbatim
  cd GeographicLib-1.30 \endverbatim
- Edit \verbatim
  include/GeographicLib/Config.h \endverbatim
  If your C++ compiler does not recognize the long double type
  (unlikely), insert \code
  #undef HAVE_LONG_DOUBLE \endcode
  If your machine using big endian ordering, then insert \code
  #define WORDS_BIGENDIAN 1 \endcode
- Build and install the software: \verbatim
  make             # compile the library and the examples
  make install     # as root \endverbatim
  The installation is in directories under /usr/local.  You
  can specify a different installation directory with, for example,
  \verbatim
  make PREFIX=/tmp/geographic install \endverbatim
- The headers, library, and utilities are installed in the
  include/GeographicLib, lib, and bin directories under
  <code>PREFIX</code>.  This installation method does not compile
  the Matlab/Octave interface;  however the source for the interface is
  installed in libexec/GeographicLib/matlab, see \ref matlab of
  instructions on compiling the interface.  For documentation, open
  <a href="index.html">
  share/doc/GeographicLib/html/index.html</a> in a web browser.

\section windows Installation on Windows

This method requires Visual Studio 2005, 2008, or 2010.  This builds a
static library and the utilities.  If you only have Visual Studio 2003,
use cmake to create the necessary solution file, see \ref cmake.  (cmake
is needed to build the Matlab interface and to run the tests.)
- Unpack the source, running \verbatim
  unzip -q GeographicLib-1.30.zip \endverbatim
- Open GeographicLib-1.30/windows/GeographicLib-vc10.sln in Visual Studio
  2010 (for Visual Studio 2005 and 2008, replace -vc10 by -vc8 or -vc9).
- Pick the build type (e.g., Release), and select "Build Solution".
- The library and the compiled examples are in the windows/Release.
- Copy the library windows/Release/Geographic.lib and the
  headers in include/GeographicLib somewhere convenient.  The
  headers should remain in a directory named %GeographicLib.  If you
  expect to use the Matlab/Octave interface, copy matlab/*.m and
  matlab/*.cpp to a directory in your matlab/octave path, see \ref
  matlab for instructions on compiling the interface.  For documentation,
  open
  <a href="index.html">doc/html/index.html</a> in a web
  browser.

\section windowsbin Using a binary installer for Windows

Use this method if you only need to use the %GeographicLib utilities.
The header files and static library are provided, but can only be used
by Visual Studio 2010 in release mode.   However, if you plan to use the
library, it is advisable to build it with the compiler you are using for
your own code using either \ref cmake or \ref windows.

Download and run
<a href="http://sourceforge.net/projects/geographiclib/files/distrib/GeographicLib-1.30-win32.exe">
GeographicLib-1.30-win32.exe</a> or
<a href="http://sourceforge.net/projects/geographiclib/files/distrib/GeographicLib-1.30-win64.exe">
GeographicLib-1.30-win64.exe</a>:
 - read the MIT/X11 License agreement,
 - select whether you want your PATH modified,
 - select the installation folder, by default
   C:\\pkg-vc10\\GeographicLib-1.30 or C:\\pkg-vc10-x64\\GeographicLib-1.30,
 - select the start menu folder,
 - and install.
 .
(Note that the default installation folder adheres the the convention
given in \ref cmake.)  The start menu will now include links to the
documentation for the library and for the utilities (and a link for
uninstalling the library).  If you ask for your PATH to be modified, it
will include C:/pkg-vc10/GeographicLib-1.30/bin where the utilities are
installed.  The headers and library are installed in the
include/GeographicLib and lib folders.  With the 64-bit installer, the
Matlab interface is installed in the matlab folder.  Add this to your
path in Matlab to access this interface.  The binaries were built
using Visual Studio 10 (2010) in Release mode.  The Matlab interface
was compiled with Matlab R2010b 64-bit, however it may work with some
other 64-bit versions of Matlab.

\section qt Building the library for use with Qt

If Qt is using a standard compiler, then build %GeographicLib with that
same compiler (and optimization flags) as Qt.

If you are using the mingw compiler on Windows for Qt, then you need to
build %GeographicLib with mingw.  You can accomplish this with cmake
under cygwin with, for example, \verbatim
  export PATH="`cygpath -m c:/QtSDK/mingw/bin`:$PATH"
  mkdir BUILD
  cd BUILD
  cmake -G "MinGW Makefiles" -D CMAKE_INSTALL_PREFIX=C:/pkg-mingw/GeographicLib ..
  mingw32-make
  mingw32-make install \endverbatim
If cmake complains that sh mustn't be in your path, invoke cmake with
\verbatim
  env PATH="$( echo $PATH | tr : '\n' |
  while read d; do test -f "$d/sh.exe" || echo -n "$d:"; done |
  sed 's/:$//' )" \
  cmake -G "MinGW Makefiles" -D CMAKE_INSTALL_PREFIX=C:/pkg-mingw/GeographicLib ..
\endverbatim
If cmake is not available, there is a simple project file that compiles
the %GeographicLib library <i>only</i> with the Qt compiler: \verbatim
  cd src
  qmake Geographic.pro   # configure the project
  make                   # build the library \endverbatim
The library will be in the src directory (or the src/release or
src/debug directory under Windows).  The headers are in
include/GeographicLib.

\section maintainer Maintainer tasks

Check the code out of git with \verbatim
  git clone -b master git://git.code.sf.net/p/geographiclib/code geographiclib
\endverbatim
Here the "master" branch is checked out.  There are three branches in
the git repository:
- <b>master</b>: the main branch for code maintainence.  Releases are
  tagged on this branch as, e.g., v1.30.
- <b>devel</b>: the development branch; changes made here are merged
  into master.
- <b>release</b>: the release branch created by unpacking the source
  releases (git is \e not used to merge changes from the other
  branches into this branch).  This is the \e default branch of the
  repository (the branch you get if cloning the repository without
  specifying a branch).  This differs from the master branch in that
  some administrative files are excluded while some intermediate files
  are included (in order to aid building on as many platforms as
  possible).  Releases are tagged on this branch as, e.g., r1.30.
.
The autoconf configuration script and the formatted man pages are not
checked into master branch of the repository.  In order to create the
autoconf configuration script, run \verbatim
  sh autogen.sh \endverbatim
in the top level directory.  Provided you are running on a system with
doxygen, pod2man, and pod2html installed, then you can create the
documentation and the man pages by building the system using cmake or
configure.

In the case of cmake, you then run \verbatim
  make dist \endverbatim
which will copy the man pages from the build directory back into the
source tree and package the resulting source tree for distribution as
\verbatim
  GeographicLib-1.30.tar.gz
  GeographicLib-1.30.zip \endverbatim
Finally, \verbatim
  make package \endverbatim
or building PACKAGE in Visual Studio will create a binary installer for
%GeographicLib.  For Windows, this requires in the installation of
<a href="http://nsis.sourceforge.net">NSIS</a>.

With configure, run \verbatim
  make dist-gzip \endverbatim
which will create the additional files and packages the results ready
for distribution as \verbatim
  geographiclib-1.30.tar.gz \endverbatim

<center>
Back to \ref intro.  Forward to \ref start.  Up to \ref contents.
</center>
**********************************************************************/
/**
\page start Getting started

<center>
Back to \ref install.  Forward to \ref utilities.  Up to \ref contents.
</center>

Much (but not all!) of the useful functionality of %GeographicLib is
available via simple command line utilities.  Interfaces to some of them
are available via the web.  See \ref utilities for documentation on
these.

In order to use %GeographicLib from C++ code, you will need to
- Include the header files for the %GeographicLib classes in your code.
  E.g., \code
  #include <GeographicLib/LambertConformalConic.hpp> \endcode
- Include the GeographicLib:: namespace prefix to the %GeographicLib classes,
  or include \code
  using namespace GeographicLib; \endcode
  in your code.
- Tell the compiler where to find the header files.  With g++ and with
  /usr/local specified as the installation directory,
  this is accomplished with \verbatim
  g++ -c -g -O3 -I/usr/local/include testprogram.cpp
  \endverbatim
  With Visual Studio, specify the include directory in the IDE via,
  e.g.,
  \verbatim
  C/C++ -> General -> Additional Include Directories = C:\pkg-vc10\GeographicLib\include
  \endverbatim
- Tell the linker the name, Geographic, and location of the
  library.  Using g++ as the linker, you would use \verbatim
  g++ -g -o testprogram testprogram.o -L/usr/local/lib -lGeographic
  \endverbatim
  With Visual Studio, you supply this information in the IDE via,
  e.g.,
  \verbatim
  Linker -> Input -> Additional Dependencies = Geographic.lib
  Linker -> General -> Additional Library Directories = C:\pkg-vc10\GeographicLib\lib
  \endverbatim
  Note that the library name is <b>Geographic</b> and not
  %GeographicLib.  If the library was configured and built with cmake,
  then the debug version of the library is called Geographic_d.lib.
- Tell the runtime environment where to find the shared library
  (assuming you compiled %GeographicLib as a shared library).  With g++,
  this is accomplished by modifying the link line above to read \verbatim
  g++ -g -o testprogram testprogram.o -Wl,-rpath=/usr/local/lib -L/usr/local/lib -lGeographic
  \endverbatim
  (There are two other ways to specify the location of shared libraries
  at runtime: (1) define the environment variable
  <code>LD_LIBRARY_PATH</code> to be a colon-separated list of
  directories to search; (2) as <b>root</b>, specify /usr/local/lib as a
  directory searched by ldconfig(8).)  On Windows, you need to ensure
  that Geographic.dll is in the same directory as your executable or
  else include the directory containing the dll in your
  <code>PATH</code>.
- If you're using cmake to configure and build your project, then
  instead of the previous three steps, insert \verbatim
  # Put in your top-level CMakeLists.txt
  find_package (GeographicLib 1.9 REQUIRED)
  # In the top-level CMakeLists.txt or wherever your code is compiled
  include_directories (${GeographicLib_INCLUDE_DIRS})
  # In the CMakeLists.txt where you define your executable targets
  target_link_libraries (program1 ${GeographicLib_LIBRARIES})
  target_link_libraries (program2 ${GeographicLib_LIBRARIES}) \endverbatim
  in your CMakeLists.txt files (as noted).  find_package command should
  find the library and set the required cmake variables (including
  <code>GeographicLib_FOUND = 1</code>).  For find_package to be able
  to locate %GeographicLib, it may be necessary for %GeographicLib to be
  built and installed with cmake (see \ref cmake) instead of using some
  other installation method.  In addition, on Windows systems, you
  should specify <code>CMAKE_PREFIX_PATH</code> with, for example,
  \verbatim
  cmake -G "Visual Studio 10" -D CMAKE_PREFIX_PATH=C:/pkg-vc10 -D CMAKE_INSTALL_PREFIX=C:/pkg-vc10/XYZProject ..
  \endverbatim
  You can direct find_package to a specific version of %GeographicLib
  by specifying the cmake variable <code>GeographicLib_DIR</code> to
  point to the directory containing the configuration file
  <code>geographiclib-config.cmake</code>; this should normally the
  cmake directory in the installation tree, but can also be a build
  directory for %GeographicLib.  On Windows systems, find_package will
  only look for versions of %GeographicLib built with the right
  version of the compiler.  (If you used a non-cmake method of
  installing %GeographicLib, you can try copying
  cmake/FindGeographicLib.cmake to somewhere in your
  <code>CMAKE_MODULE_PATH</code> in order for find_package to work.
  However, this method has not been thoroughly tested.)
- For Windows, ensure that you build the library with the same version
  of the compiler and you use for your code.  In addition, you have
  to build a separate debug version of the library when you can
  compiling your code in debug mode.  The configuration of cmake
  described above ensures that the right version of the compiler is used
  and, assuming you built and installed %GeographicLib in debug mode, the
  debug version of the library Geographic_d.lib will be used.

Here is a very simple test code, which uses the GeographicLib::Geodesic
class:
\include example-Geodesic-small.cpp
This example is <code>examples/example-Geodesic-small.cpp</code>.  If you
compile, link, and run it according to the instructions above, it should
print out \verbatim
  5551.76 km
\endverbatim

The next steps are:
 - Learn about and run the \ref utilities.
 - Read the section, \ref organization, for an overview of the library.
 - Browse the <a href="annotated.html">Class List</a> for full documentation
   on the classes in the library.
 - Look at the example code in the examples directory.  Each file
   provides a very simple standalone example of using one %GeographicLib
   class.  These are included in the descriptions of the classes.
 - Look at the source code for the utilities in the tools directory for
   more examples of using %GeographicLib from C++ code, e.g.,
   GeodesicProj.cpp is a program to performing various geodesic
   projections.

Here's a list of some of the abbreviations used here with links to the
corresponding Wikipedia articles:
 - <a href="http://en.wikipedia.org/wiki/WGS84">
   WGS84</a>, World Geodetic System 1984.
 - <a href="http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system">
   UTM</a>, Universal Transverse Mercator coordinate system.
 - <a href="http://en.wikipedia.org/wiki/Universal_Polar_Stereographic">
   UPS</a>, Universal Polar Stereographic coordinate system.
 - <a href="http://en.wikipedia.org/wiki/Military_grid_reference_system">
   MGRS</a>, Military Grid Reference System.
 - <a href="http://en.wikipedia.org/wiki/Geoid">
   EGM</a>, Earth Gravity Model.
 - <a href="http://en.wikipedia.org/wiki/World_Magnetic_Model">
   WMM</a>, World Magnetic Model.
 - <a href="http://en.wikipedia.org/wiki/IGRF">
   IGRF</a>, International Geomagnetic Reference Field.

<center>
Back to \ref install.  Forward to \ref utilities.  Up to \ref contents.
</center>
**********************************************************************/
/**
\page utilities Utility programs

<center>
Back to \ref start.  Forward to \ref organization.  Up to \ref contents.
</center>

Various utility programs are provided with %GeographicLib.  These should
be installed in a directory included in your PATH (e.g.,
/usr/local/bin).  These programs are wrapper programs that invoke
the underlying functionality provided by the library.

The utilities are
 - <a href="GeoConvert.1.html">
   <b>GeoConvert</b></a>: convert geographic coordinates using
   GeographicLib::GeoCoords.
 - <a href="GeodSolve.1.html">
   <b>GeodSolve</b></a>: perform geodesic calculations using
   GeographicLib::Geodesic and GeographicLib::GeodesicLine.
 - <a href="Planimeter.1.html">
   <b>Planimeter</b></a>: compute the area of geodesic polygons using
   GeographicLib::PolygonArea.
 - <a href="TransverseMercatorProj.1.html">
   <b>TransverseMercatorProj</b></a>: convert between geographic
   and transverse Mercator.  This is for testing
   GeographicLib::TransverseMercatorExact and
   GeographicLib::TransverseMercator.
 - <a href="CartConvert.1.html">
   <b>CartConvert</b></a>: convert geodetic coordinates to geocentric or
   local cartesian using GeographicLib::Geocentric and
   GeographicLib::LocalCartesian.
 - <a href="GeodesicProj.1.html">
   <b>GeodesicProj</b></a>: perform projections based on geodesics
   using GeographicLib::AzimuthalEquidistant, GeographicLib::Gnomonic,
   and GeographicLib::CassiniSoldner.
 - <a href="ConicProj.1.html">
   <b>ConicProj</b></a>: perform conic projections using
   GeographicLib::LambertConformalConic and
   GeographicLib::AlbersEqualArea.
 - <a href="GeoidEval.1.html">
   <b>GeoidEval</b></a>: look up geoid heights using
   GeographicLib::Geoid.
 - <a href="Gravity.1.html">
   <b>Gravity</b></a>: compute the earth's gravitational field using
   GeographicLib::GravityModel and GeographicLib::GravityCircle.
 - <a href="MagneticField.1.html">
   <b>MagneticField</b></a>: compute the earth's magnetic field using
   GeographicLib::MagneticModel and GeographicLib::MagneticCircle.
 .
The documentation for these utilities is in the form of man pages.  This
documentation can be accessed by clicking on the utility name in the
list above, running the man command on Unix-like systems, or by invoking
the utility with the --help option.  A brief summary of usage is given
by invoking the utility with the -h option.  The version of the utility
is given by the --version option.

The utilities all accept data on standard input, transform it in some
way, and print the results on standard output.  This makes the utilities
easy to use within scripts to transform tabular data; however they can
also be used interactively, often with the input supplied via a pipe,
e.g.,
 - echo 38SMB4488 | GeoConvert -d

Online versions of four of these utilities are provided:
 - <a href="http://geographiclib.sf.net/cgi-bin/GeoConvert">GeoConvert</a>
 - <a href="http://geographiclib.sf.net/cgi-bin/GeodSolve">GeodSolve</a>
 - <a href="http://geographiclib.sf.net/cgi-bin/Planimeter">Planimeter</a>
 - <a href="http://geographiclib.sf.net/cgi-bin/GeoidEval">GeoidEval</a>

<center>
Back to \ref start.  Forward to \ref organization.  Up to \ref contents.
</center>
**********************************************************************/
/**
\page organization Code organization

<center>
Back to \ref utilities.  Forward to \ref other.  Up to \ref contents.
</center>

Here is a brief description of the relationship between the various
components of %GeographicLib.  All of these are defined in the
GeographicLib namespace.

GeographicLib::TransverseMercator, GeographicLib::PolarStereographic,
GeographicLib::LambertConformalConic, and GeographicLib::AlbersEqualArea
provide the basic projections.  The constructors for these classes
specify the ellipsoid and the forward and reverse projections are
implemented as const member functions.  TransverseMercator uses
Kr&uuml;ger's series which have been extended to sixth order in the
square of the eccentricity.  PolarStereographic, LambertConformalConic,
and AlbersEqualArea use the exact formulas for the projections (e.g.,
from Snyder).

GeographicLib::TransverseMercator::UTM and
GeographicLib::PolarStereographic::UPS are const static instantiations
specific for the WGS84 ellipsoid with the UTM and UPS scale factors.
(These do \e not add the standard false eastings or false northings for
UTM and UPS.)  Similarly GeographicLib::LambertConformalConic::Mercator
is a const static instantiation of this projection for a WGS84 ellipsoid
and a standard parallel of 0 (which gives the Mercator projection).
GeographicLib::AlbersEqualArea::CylindricalEqualArea,
AzimuthalEqualAreaNorth, and AzimuthalEqualAreaSouth, likewise provide
special cases of the equal area projection.

GeographicLib::UTMUPS uses TransverseMercator::UTM and
PolarStereographic::UPS to perform the UTM and UPS
projections.  The class offers a uniform interface to UTM and UPS by
treating UPS as UTM zone 0.  This class stores no internal state and the
forward and reverse projections are provided via static member
functions.  The forward projection offers the ability to override the
standard UTM/UPS choice and the UTM zone.

GeographicLib::MGRS transforms between UTM/UPS coordinates and MGRS.
UPS coordinates are handled as UTM zone 0.  This class stores no
internal state and the forward (UTM/UPS to MGRS) and reverse (MGRS to
UTM/UPS) conversions are provided via static member functions.

GeographicLib::GeoCoords holds a single geographic location which may be
specified as latitude and longitude, UTM or UPS, or MGRS.  Member
functions are provided to convert between coordinate systems and to
provide formatted representations of them.
<a href="GeoConvert.1.html">GeoConvert</a> is a simple command line
utility to provide access to the GeoCoords class.

GeographicLib::TransverseMercatorExact is a drop in replacement for
TransverseMercator which uses the exact formulas, based on elliptic
functions, for the projection as given by Lee.
<a href="TransverseMercatorProj.1.html">TransverseMercatorProj</a> is a
simple command line utility to test to the TransverseMercator and
TransverseMercatorExact.

GeographicLib::Geodesic and GeographicLib::GeodesicLine perform geodesic
calculations.  The constructor for GeographicLib::Geodesic specifies the
ellipsoid and the direct and inverse calculations are implemented as
const member functions.  GeographicLib::Geocentric::WGS84 is a const
static instantiation of Geodesic specific for the WGS84 ellipsoid.  In
order to perform a series of direct geodesic calculations on a single
line, the GeographicLib::GeodesicLine class can be used.  This packages
all the information needed to specify a geodesic.  A const member
function returns the coordinates a specified distance from the starting
point.  <a href="GeodSolve.1.html">GeodSolve</a> is a simple command
line utility to perform geodesic calculations.
GeographicLib::PolygonArea is a class which compute the area of geodesic
polygons using the Geodesic class and <a
href="Planimeter.1.html">Planimeter</a> is a command line utility for
the same purpose.  GeographicLib::AzimuthalEquidistant,
GeographicLib::CassiniSoldner, and GeographicLib::Gnomonic are
projections based on the Geodesic class.  <a
href="GeodesicProj.1.html">GeodesicProj</a> is a command line utility to
exercise these projections.

GeographicLib::GeodesicExact and GeographicLib::GeodesicLineExact are
drop in replacements for GeographicLib::Geodesic and
GeographicLib::GeodesicLine in which the solution is given in terms of
elliptic integrals (computed by GeographicLib::EllipticFunction).
These classes should be used if the absolute value of the flattening
exceeds 0.02.  The -E option to <a href="GeodSolve.1.html">GeodSolve</a>
uses these classes.

GeographicLib::Geocentric and GeographicLib::LocalCartesian convert between
geodetic and geocentric or a local cartesian system.  The constructor for
specifies the ellipsoid and the forward and reverse projections are
implemented as const member functions.  GeographicLib::Geocentric::WGS84 is a
const static instantiation of Geocentric specific for the WGS84 ellipsoid.
<a href="CartConvert.1.html">CartConvert</a> is a simple command line
utility to provide access to these classes.

GeographicLib::Geoid evaluates geoid heights by interpolation.  This is
provided by the operator() member function.
<a href="GeoidEval.1.html">GeoidEval</a> is a simple command line
utility to provide access to this class.  This class requires
installation of data files for the various geoid models; see \ref
geoidinst for details.

GeographicLib::Ellipsoid is a class which performs latitude
conversions and returns various properties of the ellipsoid.

GeographicLib::GravityModel evaluates the earth's gravitational field
using a particular gravity model.  Various member functions return the
gravitational field, the gravity disturbance, the gravity anomaly, and
the geoid height <a href="Gravity.1.html">Gravity</a> is a simple
command line utility to provide access to this class.  If the field
several points on a circle of latitude are sought then use
GeographicLib::GravityModel::Circle to return a
GeographicLib::GravityCircle object whose member functions performs the
calculations efficiently.  (This is particularly important for high
degree models such as EGM2008.)  These classes requires installation of
data files for the various gravity models; see \ref gravityinst for
details.

GeographicLib::MagneticModel evaluates the earth's magnetic field using
a particular magnetic model.  The field is provided by the operator()
member function.  <a href="MagneticField.1.html">MagneticField</a> is a
simple command line utility to provide access to this class.  If the
field several points on a circle of latitude are sought then use
GeographicLib::MagneticModel::Circle to return a
GeographicLib::MagneticCircle object whose operator() member function
performs the calculation efficiently.  (This is particularly important
for high degree models such as emm2010.)  These classes requires
installation of data files for the various magnetic models; see \ref
magneticinst for details.

GeographicLib::Constants, GeographicLib::Math, GeographicLib::Utility,
GeographicLib::DMS, are general utility class which are used
internally by the library; in addition GeographicLib::EllipticFunction
is used by GeographicLib::TransverseMercatorExact,
GeographicLib::Ellipsoid, and GeographicLib::GeodesicExact, and
GeographicLib::GeodesicLineExact; GeographicLib::Accumulator is used
by GeographicLib::PolygonArea, and GeographicLib::SphericalEngine,
GeographicLib::CircularEngine, GeographicLib::SphericalHarmonic,
GeographicLib::SphericalHarmonic1, and
GeographicLib::SphericalHarmonic2 facilitate the summation of
spherical harmonic series which is needed by and
GeographicLib::MagneticModel and GeographicLib::MagneticCircle.  One
important definition is GeographicLib::Math::real which is the type
used for real numbers.  This allows the library to be easily switched
to using floats, doubles, or long doubles.  However all the testing
has been with real set to double and the library should be installed
in this way.

In general, the constructors for the classes in %GeographicLib check
their arguments and throw GeographicLib::GeographicErr exceptions with a
explanatory message if these are illegal.  The member functions, e.g.,
the projections implemented by TransverseMercator and
PolarStereographic, the solutions to the geodesic problem, etc.,
typically do <i>not</i> check their arguments; the calling program
should ensure that the arguments are legitimate.  However, the functions
implemented by UTMUPS, MGRS, and GeoCoords do check their arguments and
may throw GeographicLib::GeographicErr exceptions.  Similarly Geoid may
throw exceptions on file errors.  If a function does throw an exception,
then the function arguments used for return values will not have been
altered.

%GeographicLib attempts to act sensibly with NaNs.  NaNs in constructors
typically throw errors (an exception is GeodesicLine).  However, calling
the class functions with NaNs as arguments is not an error; NaNs are
returned as appropriate.  "INV" is treated as an invalid zone
designation by UTMUPS.  "INVALID" is the corresponding invalid MGRS
string.  (Similary "nan" is the equivalent invalid Geohash.)  NaNs allow
the projection of polylines which are separated by NaNs; in this format
they can be easily plotted in Matlab.

A note about portability.  For the most part, the code uses standard C++
and should be able to be deployed on any system with a modern C++
compiler.  System dependencies come into
 - GeographicLib::Math -- GeographicLib needs to define functions such
   as atanh for systems that lack them.  The system dependence will
   disappear with the adoption of C++11 because the needed functions are
   part of that standard.
 - use of long double -- the type is used only for testing.  On systems
   which lack this data type the cmake and autoconf configuration
   methods should detect its absence and omit the code that depends on
   it.
 - GeographicLib::Geoid, GeographicLib::GravityModel, and
   GeographicLib::MagneticModel -- these class uses system-dependent
   default paths for looking up the respective datasets.  It also relies
   on getenv to find the value of the environment variables.
 - GeographicLib::Utility::readarray reads numerical data from binary
   files.  This assumes that floating point numbers are in IEEE format.
   It attempts to handled the "endianness" of the target platform using
   the WORDS_BIGENDIAN macro (which sets the compile-time constant
   GeographicLib::Math::bigendian).

<center>
Back to \ref utilities.  Forward to \ref other.  Up to \ref contents.
</center>
**********************************************************************/
/**
\page other Other interfaces (C, Fortran, Javascript, Python, Matlab, Maxima)

<center>
Back to \ref organization.  Forward to \ref geoid.  Up to \ref contents.
</center>

Interfaces to subsets of %GeographicLib are available in other languages
- \ref c-fortran
- \ref javascript.
- \ref python.
- \ref matlab.
- \ref maxima.

\section c-fortran C and Fortran

The directories <code>legacy/C</code> and <code>legacy/Fortran</code>
contain implementations of GeographicLib::Geodesic,
GeographicLib::GeodesicLine, and GeographicLib::PolygonArea in C and
Fortran respectively.  These are intended for use in old codes written
in these languages and should work any reasonably modern compiler.
These implementations are entirely self-contained and do not depend on
the rest of %GeographicLib.  Sample main programs to solve the direct
and inverse geodesic problems and to compute polygonal areas are
provided.

For documentation, see
 - <a href="C/index.html">C library for geodesics</a>,
 - <a href="Fortran/index.html">Fortran library for geodesics</a>.

\section javascript Javascript interface

The directory doc/scripts/GeographicLib contains the classes
- GeographicLib::Math
- GeographicLib::Accumulator
- GeographicLib::Geodesic
- GeographicLib::GeodesicLine
- GeographicLib::PolygonArea
- GeographicLib::DMS
.
translated into Javascript.  See Interface.js for a simple Javascript
interface to these routines (documented near the top of the file).
Examples of using this interface are
- a <a href="../scripts/geod-calc.html">geodesic calculator</a> showing
  the solution of direct and inverse geodesic problem, finding
  intermediate points on a geodesic line, and computing the area of a
  geodesic polygon.
- <a href="../scripts/geod-google.html">displaying geodesics in Google
  Maps</a> which shows the geodesic, the geodesic circle, and various
  geodesic envelopes.
.
These examples include a "stripped" version of the Javascript code, \code
<script type="text/javascript"
        src="http://geographiclib.sf.net/scripts/geographiclib.js">
</script> \endcode
which loads faster.

\section python Python interface

A python implementation of the geodesic routines from GeographicLib
are provided in the python/geographiclib directory (which is installed
as PREFIX/lib/python/site-packages/geographiclib, if
COMMON_INSTALL_PATH is ON, and as PREFIX/python/geographiclib,
otherwise).  This contains implementations of the classes
- GeographicLib::Math
- GeographicLib::Accumulator
- GeographicLib::Geodesic
- GeographicLib::GeodesicLine
- GeographicLib::PolygonArea
.
You can also download the python interface independent of the rest of
%GeographicLib from
- <a href="http://pypi.python.org/pypi/geographiclib">
  http://pypi.python.org/pypi/geographiclib</a>
.
and then unpack the .tar.gz or .zip file.

You can "install" these routines, so that they are in python's default
path with, for example \verbatim
  cd geographiclib-1.16
  python setup.py install
\endverbatim
(this will require root privileges).  Or else you can set the path
within python using \code
>>> import sys
>>> sys.path.append("/usr/local/lib/python/site-packages")
\endcode

An example of using this interface is \code
>>> from geographiclib.geodesic import Geodesic
>>> # The geodesic inverse problem
... Geodesic.WGS84.Inverse(-41.32, 174.81, 40.96, -5.50)
>>> # The geodesic direct problem
... Geodesic.WGS84.Direct(40.6, -73.8, 45, 10000e3)
>>> # How to obtain several points along a geodesic
... line = Geodesic.WGS84.Line(40.6, -73.8, 45)
>>> line.Position( 5000e3)
>>> line.Position(10000e3)
>>> # Computing the area of a geodesic polygon
... def p(lat,lon): return {'lat': lat, 'lon': lon}
...
>>> Geodesic.WGS84.Area([p(0, 0), p(0, 90), p(90, 0)])
>>> # Introductory help
... help(Geodesic)
\endcode

(Note: The initial version of setup.py was provided by Andrew MacIntyre
of the Australian Communications and Media Authority.)

\section matlab Matlab and Octave interfaces

The matlab directory contains
 - Native Matlab implementations of the geodesic routines.  To use
   these, start Matlab or Octave and run one of (for example) \verbatim
   addpath /usr/local/libexec/GeographicLib/matlab
   addpath 'C:/pkg-vc10-x64/GeographicLib-1.30/libexec/GeographicLib/matlab'
   \endverbatim
   The available functions are:
   - geoddoc: briefly descibe the routines
   - geodreckon: solve the direct geodesic problem
   - geoddistance: solve the inverse geodesic problem
   - geodarea: compute the area of ellipsoidal polygons
   .
   Use the help function to get documentation, e.g., \code
   help geoddistance \endcode
   to obtain documentation.  These functions are also available as a
   standalone package from
   <a href="http://www.mathworks.com/matlabcentral/fileexchange/">
   Matlab File Exchange</a> using the link
   - <a href="http://www.mathworks.com/matlabcentral/fileexchange/39108">
     http://www.mathworks.com/matlabcentral/fileexchange/39108</a>
 - Native Matlab implementations of projections which are related to
   geodesics.  These are
   - geodproj: briefly descibe the routines
   - eqdazim_fwd, eqdazim_inv: azimuthal equidistant
   - cassini_fwd, cassini_inv: Cassini-Soldner
   - tranmerc_fwd, tranmerc_inv: transverse Mercator
   - gnomonic_fwd, gnomonic_inv: ellipsoidal gnomonic
   - utm_fwd, utm_inv: universal transverse Mercator
   .
   These functions are also available as a package from
   <a href="http://www.mathworks.com/matlabcentral/fileexchange/">
   Matlab File Exchange</a> using the link
   - <a href="http://www.mathworks.com/matlabcentral/fileexchange/39366">
     http://www.mathworks.com/matlabcentral/fileexchange/39366</a>
   .
   (This requires that the previous package also be installed.)
 - Interface code so that some %GeographicLib classes can be accessed
   from Matlab or Octave.  The rest of this section describes how to
   compile and use these interfaces.

There are two ways of compiling the interface code: (1) using cmake and
(2) invoking the compiler from Matlab.
 - <b>Using cmake:</b> Before running cmake, configure MATLAB on
   Windows to use the same compiler that you're going to use for compiling
   %GeographicLib.  For example \verbatim
   mex.bat -setup

   Please choose your compiler for building external interface (MEX) files:
   Would you like mex to locate installed compilers [y]/n? y
   Select a compiler:
   [1] Lcc-win32 C 2.4.1 in C:\PROGRA~1\MATLAB\R2010b\sys\lcc
   [2] Microsoft Visual C++ 2010 in C:\Program Files\Microsoft Visual Studio 10.0
   [0] None

   Compiler: 2
   etc. \endverbatim
   (This will require that mex.bat is in your PATH.  With Linux, use
   <code>mex -setup</code>.)  Then configure cmake with, for example
   \verbatim
   cmake -G "Visual Studio 10" -D MATLAB_COMPILER=mex .. \endverbatim
   If you are running a 64-bit version of Matlab, be sure to select a
   64-bit generator with cmake, e.g., "Visual Studio 10 Win64".  Finally
   compile %GeographicLib with Visual Studio.  (The binary installer for
   64-bit Windows includes the compiled interface built with Visual
   Studio 10 and Matlab R2010b 64-bit).<br>
   On Linux systems, you can compile the interface for use with octave
   instead by using \verbatim
   cmake -D MATLAB_COMPILER=mkoctfile .. \endverbatim
 - <b>Invoking the compiler from Matlab or Octave:</b> Start Matlab or
   Octave and run, e.g., \code
   mex -setup
   cd 'C:/pkg-vc10-x64/GeographicLib-1.30/matlab'
   help geographiclibinterface
   geographiclibinterface('C:/pkg-vc10/GeographicLib-1.30');
   addpath(pwd);
   \endcode
   The first command allows you to select the compiler to use (which
   should be the same as that used to compile %GeographicLib).

To use the interface routines for %GeographicLib, run one of (for
example) \verbatim
  addpath /usr/local/libexec/GeographicLib/matlab
  addpath 'C:/pkg-vc10-x64/GeographicLib-1.30/libexec/GeographicLib/matlab'
\endverbatim
in Octave or Matlab.  The available functions are:
 - geodesicdirect: solve direct geodesic problem
   (see GeographicLib::Geodesic::Direct)
 - geodesicinverse: solve inverse geodesic problem
   (see GeographicLib::Geodesic::Inverse)
 - geodesicline: compute points along a geodesic
   (see GeographicLib::GeodesicLine::Position)
 - polygonarea: compute area of a geodesic polygon
   (see GeographicLib::PolygonArea)
 - utmupsforward: convert geographic coordinates to UTM/UPS
   (see GeographicLib::UTMUPS::Forward)
 - utmupsreverse: convert UTM/UPS coordinates to geographic
   (see GeographicLib::UTMUPS::Reverse)
 - mgrsforward: convert UTM/UPS coordinates to MGRS
   (see GeographicLib::MGRS::Forward)
 - mgrsreverse: convert MGRS coordinates to UTM/UPS
   (see GeographicLib::MGRS::Reverse)
 - geoidheight: compute geoid height
   (see GeographicLib::Geoid::operator()())
 - geocentricforward: convert geographic coordinates to geocentric
   (see GeographicLib::Geocentric::Forward)
 - geocentricreverse: convert geocentric coordinates to geographic
   (see GeographicLib::Geocentric::Reverse)
 - localcartesianforward: convert geographic coordinates to local cartesian
   (see GeographicLib::LocalCartesian::Forward)
 - localcartesianreverse: convert local cartesian coordinates to geographic
   (see GeographicLib::LocalCartesian::Reverse)
 .
These routines just offer a simple interface to the corresponding C++
class.  Use the help function to get documentation, e.g., \code
  help geodesicdirect \endcode
Unfortunately, the help function does not work for compiled functions in
Octave; in this case, just list the .m file, e.g., \code
  type geodesicdirect \endcode
Other useful functions, e.g., to convert from geographic
coordinates to MGRS can easily be written with Matlab code.

Note that geoidheight, when compiled with Visual Studio 2008 causes
Matlab to crash.  (The problem does not occur with Visual Studio 2005 or
Visual Studio 2010.)

\section maxima Maxima routines

Maxima is a free computer algebra system which can be downloaded from
http://maxima.sf.net.  Maxima was used to generate the series used by
GeographicLib::TransverseMercator (<a href="tmseries.mac">
tmseries.mac</a>) and GeographicLib::Geodesic (<a href="geod.mac">
geod.mac</a>) and to generate accurate data for testing
(<a href="tm.mac"> tm.mac</a> and <a href="geodesic.mac">
geodesic.mac</a>).  The latter uses Maxima's bigfloat arithmetic
together with series extended to high order or solutions in terms of
elliptic integrals (<a href="ellint.mac"> ellint.mac</a>).  These files
contain brief instructions on how to use them.

<center>
Back to \ref organization.  Forward to \ref geoid.  Up to \ref contents.
</center>
**********************************************************************/
/**
\page geoid Geoid height

<center>
Back to \ref other.  Forward to \ref gravity.  Up to \ref contents.
</center>

The gravitational equipotential surface approximately coinciding with
mean sea level is called the geoid.  The GeographicLib::Geoid class and
the <a href="GeoidEval.1.html">GeoidEval</a> utility evaluate the height
of the geoid above the WGS84 ellipsoid.  This can be used to convert
heights above mean sea level to heights above the WGS84 ellipsoid.
Because the normal to the ellipsoid differs from the normal to the geoid
(the direction of a plumb line) there is a slight ambiguity in the
measurement of heights; however for heights up to 10 km this ambiguity
is only 1 mm.

The geoid is usually approximated by an "earth gravity model" (EGM).
The models published by the NGA are:
- <b>EGM84</b>:
  http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html
- <b>EGM96</b>:
  http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html
- <b>EGM2008</b>:
  http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008
.
GeographicLib::Geoid offers a uniform way to handle all 3 geoids at a
variety of grid resolutions.  (In contrast, the software tools that NGA
offers are different for each geoid, and the interpolation programs are
different for each grid resolution.  In addition these tools are written
in Fortran with is nowadays difficult to integrate with other software.)

Unlike other components of %GeographicLib, there is a appreciable error
in the results obtained (at best, the RMS error is 1 mm).  However the
class provides methods to report the maximum and RMS errors in the
results.  The class also returns the gradient of the geoid.  This can be
used to estimate the direction of a plumb line relative to the WGS84
ellipsoid.

The GeographicLib::GravityModel class calculates geoid heights using the
underlying gravity model.  This is slower then GeographicLib::Geoid but
considerably more accurate.  This class also can accurately compute all
the components of the acceleration due to gravity (and hence the
direction of plumb line).

Go to
 - \ref geoidinst
 - \ref geoidformat
 - \ref geoidinterp
 - \ref geoidcache
 - \ref testgeoid

\section geoidinst Installing the geoid datasets

The geoid heights are computed using interpolation into a rectangular
grid.  The grids are read from data files which have been are computed
using the NGA synthesis programs in the case of the EGM84 and EGM96
models and using the NGA binary gridded data files in the case of
EGM2008.  These data files are available for download:
<center>
<table>
<caption>Available geoid data files</caption>
<tr>
 <th rowspan="2">name         <th rowspan="2">geoid    <th rowspan="2">grid
 <th rowspan="2">size\n(MB)
 <th colspan="3"><center>Download Links (size, MB)</center></th>
<tr>
 <th>tar file</th>
 <th>Windows\n installer</th>
 <th>zip file</th>
<tr>
 <td>egm84-30
 <td>EGM84
 <td><center>30'</center>
 <td><center>0.6</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm84-30.tar.bz2">
 link</a> (0.5)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm84-30.exe">
 link</a> (0.8)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm84-30.zip">
 link</a> (0.5)</center>
<tr>
 <td>egm84-15
 <td>EGM84
 <td><center>15'</center>
 <td><center>2.1</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm84-15.tar.bz2">
 link</a> (1.5)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm84-15.exe">
 link</a> (1.9)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm84-15.zip">
 link</a> (2.0)</center>
<tr>
 <td>egm96-15
 <td>EGM96
 <td><center>15'</center>
 <td><center>2.1</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm96-15.tar.bz2">
 link</a> (1.5)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm96-15.exe">
 link</a> (1.9)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm96-15.zip">
 link</a> (2.0)</center>
<tr>
 <td>egm96-5
 <td>EGM96
 <td><center>5'</center>
 <td><center> 19</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm96-5.tar.bz2">
 link</a> (11)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm96-5.exe">
 link</a> (11)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm96-5.zip">
 link</a> (17)</center>
<tr>
 <td>egm2008-5
 <td>EGM2008
 <td><center>5'</center>
 <td><center> 19</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-5.tar.bz2">
 link</a> (11)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-5.exe">
 link</a> (11)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-5.zip">
 link</a> (17)</center>
<tr>
 <td>egm2008-2_5
 <td>EGM2008
 <td><center>2.5'</center>
 <td><center> 75</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-2_5.tar.bz2">
 link</a> (35)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-2_5.exe">
 link</a> (33)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-2_5.zip">
 link</a> (65)</center>
<tr>
 <td>egm2008-1
 <td>EGM2008
 <td><center>1'</center>
 <td><center>470</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-1.tar.bz2">
 link</a> (170)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-1.exe">
 link</a> (130)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-1.zip">
 link</a> (390)</center>
</table>
</center>
The "size" column is the size of the uncompressed data and it also
gives the memory requirements for caching the entire dataset using the
GeographicLib::Geoid::CacheAll method.
At a minimum you should install egm96-5 and either egm2008-1 or
egm2008-2_5.  Many applications use the EGM96 geoid, however the use of
EGM2008 is growing. (EGM84 is rarely used now.)

For Linux and Unix systems, %GeographicLib provides a shell script
geographiclib-get-geoids (typically installed in /usr/local/sbin) which
automates the process of downloading and installing the geoid data.  For
example
\verbatim
   geographiclib-get-geoids best # to install egm84-15, egm96-5, egm2008-1
   geographiclib-get-geoids -h   # for help
\endverbatim
This script should be run as a user with write access to the
installation directory, which is typically
/usr/local/share/GeographicLib (this can be overridden with the -p
flag), and the data will then be placed in the "geoids" subdirectory.

Windows users should download and run the Windows installers.  These
will prompt for an installation directory with the default being one of
\verbatim
   C:/Documents and Settings/All Users/Application Data/GeographicLib
   C:/ProgramData/GeographicLib
\endverbatim
(which you probably should not change) and the data is installed in the
"geoids" sub-directory.  (The second directory name is an alternate name
that Windows 7 for the "Application Data" directory.)

Otherwise download \e either the tar.bz2 file \e or the zip file (they
have the same contents).  If possible use the tar.bz2 files, since these
are compressed about 2 times better than the zip file.  To unpack these,
run, for example
\verbatim
   mkdir -p /usr/local/share/GeographicLib
   tar xofjC egm96-5.tar.bz2 /usr/local/share/GeographicLib
   tar xofjC egm2008-2_5.tar.bz2 /usr/local/share/GeographicLib
   etc.
\endverbatim
and, again, the data will be placed in the "geoids" subdirectory.

However you install the geoid data, all the datasets should
be installed in the same directory.  GeographicLib::Geoid and
<a href="GeoidEval.1.html">GeoidEval</a> uses a compile time default to
locate the datasets.  This is
- /usr/local/share/GeographicLib/geoids, for non-Windows systems
- C:/Documents and Settings/All Users/Application Data/GeographicLib/geoids,
  for Windows systems
.
consistent with the examples above.  This may be overridden at run-time
by defining the GEOID_PATH or the GEOGRAPHIC_DATA environment variables;
see GeographicLib::Geoid::DefaultGeoidPath() for details.  Finally, the
path may be set using the optional second argument to the
GeographicLib::Geoid constructor or with the "-d" flag to
<a href="GeoidEval.1.html">GeoidEval</a>.  Supplying the "-h" flag to
<a href="GeoidEval.1.html">GeoidEval</a> reports the default geoid path
for that utility.  The "-v" flag causes GeoidEval to report the full
path name of the data file it uses.

\section geoidformat The format of the geoid data files

The gridded data used by the GeographicLib::Geoid class is stored in
16-bit PGM files.  Thus the data for egm96-5 might be stored in the file
- /usr/local/share/GeographicLib/geoids/egm96-5.pgm
.
PGM simple graphic format with the following properties
- it is well documented
  <a href="http://netpbm.sf.net/doc/pgm.html">here</a>;
- there are no separate "little-endian" and "big-endian" versions;
- it uses 1 or 2 bytes per pixel;
- pixels can be randomly accessed;
- comments can be added to the file header;
- it is sufficiently simple that it can be easily read without using the
  <a href="http://netpbm.sf.net/doc/libnetpbm.html">libnetpbm</a>
  library (and thus we avoid adding a software dependency to
  %GeographicLib).
.
The major drawback of this format is that since there are only 65535
possible pixel values, the height must be quantized to 3 mm.  However,
the resulting quantization error (up to 1.5 mm) is typically smaller
than the linear interpolation errors.  The comments in the header for
egm96-5 are
\verbatim
   # Geoid file in PGM format for the GeographicLib::Geoid class
   # Description WGS84 EGM96, 5-minute grid
   # URL http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html
   # DateTime 2009-08-29 18:45:03
   # MaxBilinearError 0.140
   # RMSBilinearError 0.005
   # MaxCubicError 0.003
   # RMSCubicError 0.001
   # Offset -108
   # Scale 0.003
   # Origin 90N 0E
   # AREA_OR_POINT Point
   # Vertical_Datum WGS84
\endverbatim
Of these lines, the Scale and Offset lines are required and define the
conversion from pixel value to height (in meters) using \e height =
\e offset + \e scale \e pixel.  The Geoid constructor also reads the
Description, DateTime, and error lines (if present) and stores the
resulting data so that it can be returned by
GeographicLib::Geoid::Description, GeographicLib::Geoid::DateTime,
GeographicLib::Geoid::MaxError, and GeographicLib::Geoid::RMSError
methods.  The other lines serve as additional documentation but are not
used by this class.  Accompanying egm96-5.pgm (and similarly with the
other geoid data files) are two files egm96-5.wld and
egm96-5.pgm.aux.xml.  The first is an ESRI "world" file and the second
supplies complete projection metadata for use by
<a href="http://www.gdal.org">GDAL</a>.  Neither of these files is read
by GeographicLib::Geoid.

You can use gdal_translate to convert the data files to a standard
GeoTiff, e.g., with
\verbatim
   gdal_translate -ot Float32 -scale 0 65000 -108 87 egm96-5.pgm egm96-5.tif
\endverbatim
The arguments to -scale here are specific to the Offset and Scale
parameters used in the pgm file (note 65000 * 0.003 - 108 = 87).  You
can check these by running <a href="GeoidEval.1.html">GeoidEval</a> with
the "-v" option.

Here is a sample script which uses GDAL to create a 1-degree
squared grid of geoid heights at 3&quot; resolution (matching DTED1) by
bilinear interpolation.
\verbatim
   #! /bin/sh
   lat=37
   lon=067
   res=3                           # resolution in seconds
   TEMP=`mktemp junkXXXXXXXXXX`    # temporary file for GDAL
   gdalwarp -q -te `echo $lon $lat $res | awk '{
       lon = $1; lat = $2; res = $3;
       printf "%.14f %.14f %.14f %.14f",
           lon  -0.5*res/3600, lat  -0.5*res/3600,
           lon+1+0.5*res/3600, lat+1+0.5*res/3600;
   }'` -ts $((3600/res+1)) $((3600/res+1)) -r bilinear egm96-5.tif $TEMP
   gdal_translate -quiet \
       -mo AREA_OR_POINT=Point \
       -mo Description="WGS84 EGM96, $res-second grid" \
       -mo Vertical_Datum=WGS84 \
       -mo Tie_Point_Location=pixel_corner \
       $TEMP e$lon-n$lat.tif
   rm -f $TEMP
\endverbatim

Because the pgm files are uncompressed, they can take up a lot of disk
space.  Some compressed formats compress in tiles and so might be
compatible with the requirement that the data can be randomly accessed.
In particular gdal_translate can be used to convert the pgm files to
compressed tiff files with
\verbatim
gdal_translate -co COMPRESS=LZW -co PREDICTOR=2 \
   -co TILED=YES -co BLOCKXSIZE=256 -co BLOCKYSIZE=256 \
   egmyy-g.pgm egmyy-g.tif
\endverbatim
The resulting files sizes are
\verbatim
                  pgm      tif
    egm84-30      0.6 MB   0.5 MB
    egm84-15      2.1 MB   1.4 MB
    egm96-15      2.1 MB   1.5 MB
    egm96-5        19 MB   8.5 MB
    egm2008-5      19 MB   9.8 MB
    egm2008-2_5    75 MB    28 MB
    egm2008-1     470 MB    97 MB
\endverbatim
Currently, there are no plans for %GeographicLib to support this
compressed format.

\section geoidinterp Interpolating the geoid data

GeographicLib::Geoid evaluates the geoid height using bilinear or cubic
interpolation.  The gradient of the geoid height is obtained by
differentiating the interpolated height and referencing the result to
distance on the WGS84 ellipsoid.

<b>WARNING</b>: Although GeographicLib::Geoid computes the gradient of
the geoid height, this is subject to large quantization errors as a
result of the way that the geoid data is stored.  This is particularly
acute for fine grids, at high latitudes, and for the easterly gradient.
If you need to compute the direction of the acceleration due to gravity
accurately, you should use GeographicLib::GravityModel::Gravity.

The bilinear interpolation is based on the values at the 4 corners of
the enclosing cell.  The interpolated height is a continuous function of
position; however the gradient has discontinuities are cell boundaries.
The quantization of the data files exacerbates the errors in the
gradients.

The cubic interpolation is a least-squares fit to the values on a
12-point stencil with weights as follows:
\verbatim
   . 1 1 .
   1 2 2 1
   1 2 2 1
   . 1 1 .
\endverbatim
The cubic is constrained to be independent of longitude when evaluating
the height at one of the poles.  Cubic interpolation is considerably
more accurate than bilinear interpolation; however, in this
implementation there are small discontinuities in the heights are cell
boundaries.  The over-constrained cubic fit slightly reduces the
quantization errors on average.

The algorithm for the least squares fit is taken from, F. H. Lesh,
Multi-dimensional least-squares polynomial curve fitting, CACM 2, 29-30
(1959).  This algorithm is not part of GeographicLib::Geoid; instead it is
implemented as
<a href="http://en.wikipedia.org/wiki/Maxima_(software)">Maxima</a>
code which is used to precompute the matrices to convert the function
values on the stencil into the coefficients from the cubic polynomial.
This code is included as a comment in Geoid.cpp.

The interpolation methods are quick and give good accuracy.  Here is a
summary of the combined quantization and interpolation errors for the
heights.
<center><table>
<caption>Interpolation and quantization errors for geoid heights</caption>
<tr>
 <th rowspan="2">name
 <th rowspan="2">geoid
 <th rowspan="2">grid
 <th colspan="2"><center>bilinear error</center></th>
 <th colspan="2"><center>cubic error</center></th>
<tr>
 <th><center>max</center><th><center>rms</center>
 <th><center>max</center><th><center>rms</center>
<tr>
 <td>egm84-30   <td>EGM84  <td><center>30'  </center>
 <td><center>1.546 m</center><td><center>70 mm</center>
 <td><center>0.274 m</center><td><center>14 mm</center>
<tr>
 <td>egm84-15   <td>EGM84  <td><center>15'  </center>
 <td><center>0.413 m</center><td><center>18 mm</center>
 <td><center>0.021 m</center><td><center>1.2 mm</center>
<tr>
 <td>egm96-15   <td>EGM96  <td><center>15'  </center>
 <td><center>1.152 m</center><td><center>40 mm</center>
 <td><center>0.169 m</center><td><center>7.0 mm</center>
<tr>
 <td>egm96-5    <td>EGM96  <td><center> 5'  </center>
 <td><center>0.140 m</center><td><center>4.6 mm</center>
 <td><center>0.0032 m</center><td><center>0.7 mm</center>
<tr>
 <td>egm2008-5  <td>EGM2008<td><center> 5'  </center>
 <td><center>0.478 m</center><td><center>12 mm</center>
 <td><center>0.294 m</center><td><center>4.5 mm</center>
<tr>
 <td>egm2008-2_5<td>EGM2008<td><center> 2.5'</center>
 <td><center>0.135 m</center><td><center>3.2 mm</center>
 <td><center>0.031 m</center><td><center>0.8 mm</center>
<tr>
 <td>egm2008-1  <td>EGM2008<td><center> 1'  </center>
 <td><center>0.025 m</center><td><center>0.8 mm</center>
 <td><center>0.0022 m</center><td><center>0.7 mm</center>
</table></center>
The errors are with respect to the specific NGA earth gravity model
(not to any "real" geoid).  The RMS values are obtained using 5 million
uniformly distributed random points.  The maximum values are obtained by
evaluating the errors using a different grid with points at the
centers of the original grid.  (The RMS difference between EGM96 and
EGM2008 is about 0.5 m.  The RMS difference between EGM84 and EGM96 is
about 1.5 m.)

The errors in the table above include the quantization errors that arise
because the height data that GeographicLib::Geoid uses are quantized to
3 mm.  If, instead, GeographicLib::Geoid were to use data files without
such quantization artifacts, the overall error would be reduced <i>but
only modestly</i> as shown in the following table, where only the
changed rows are included and where the changed entries are given in
bold:
<center><table>
<caption>Interpolation (only!) errors for geoid heights</caption>
<tr>
 <th rowspan="2">name
 <th rowspan="2">geoid
 <th rowspan="2">grid
 <th colspan="2"><center>bilinear error</center></th>
 <th colspan="2"><center>cubic error</center></th>
<tr>
 <th><center>max</center><th><center>rms</center>
 <th><center>max</center><th><center>rms</center>
<tr>
 <td>egm96-5    <td>EGM96  <td><center> 5'  </center>
 <td><center>0.140 m</center><td><center>4.6 mm</center>
 <td><center><b>0.0026 m</b></center><td><center><b>0.1 mm</b></center>
<tr>
 <td>egm2008-2_5<td>EGM2008<td><center> 2.5'</center>
 <td><center>0.135 m</center><td><center>3.2 mm</center>
 <td><center>0.031 m</center><td><center><b>0.4 mm</b></center>
<tr>
 <td>egm2008-1  <td>EGM2008<td><center> 1'  </center>
 <td><center>0.025 m</center><td><center><b>0.6 mm</b></center>
 <td><center><b>0.0010 m</b></center><td><center><b>0.011 mm</b></center>
</table></center>

\section geoidcache Caching the geoid data

A simple way of calling Geoid is as follows
\code
   #include <GeographicLib/Geoid.hpp>
   #include <iostream>
   ...
   GeographicLib::Geoid g("egm96-5");
   double lat, lon;
   while (std::cin >> lat >> lon)
      std::cout << g(lat, lon) << "\n";
   ...
\endcode

The first call to g(lat, lon) causes the data for the stencil points (4
points for bilinear interpolation and 12 for cubic interpolation) to be
read and the interpolated value returned.  A simple 0th-order caching
scheme is provided by default, so that, if the second and subsequent
points falls within the same grid cell, the data values are not reread
from the file.  This is adequate for most interactive use and also
minimizes disk accesses for the case when a continuous track is being
followed.

If a large quantity of geoid calculations are needed, the calculation
can be sped up by preloading the data for a rectangular block with
GeographicLib::Geoid::CacheArea or the entire dataset with
GeographicLib::Geoid::CacheAll.  If the requested points lie within the
cached area, the cached data values are used; otherwise the data is read
from disk as before.  Caching all the data is a reasonable choice for
the 5' grids and coarser.  Caching all the data for the 1' grid will
require 0.5 GB of RAM and should only be used on systems with sufficient
memory.

The use of caching does not affect the values returned.  Because of the
caching and the random file access, this class is \e not normally thread
safe; i.e., a single instantiation cannot be safely used by multiple
threads.  If multiple threads need to calculate geoid heights, there are
two alternatives:
 - they should all construct thread-local instantiations.
 - GeographicLib::Geoid should be constructed with \e threadsafe = true.
   This causes all the data to be read at the time of construction (and
   if this fails, an exception is thrown), the data file to be closed
   and the single-cell caching to be turned off.  The resulting object
   may then be shared safely between threads.

\section testgeoid Test data for geoids

A test set for the geoid models is available at
 - <a href="http://sf.net/projects/geographiclib/files/testdata/GeoidHeights.dat.gz">
   GeoidHeights.dat.gz</a>
 .
This is about 11 MB (compressed).  This test set consists of a set of
500000 geographic coordinates together with the corresponding geoid
heights according to various earth gravity models.

Each line of the test set gives 6 space delimited numbers
 - latitude (degrees, exact)
 - longitude (degrees, exact)
 - EGM84 height (meters, accurate to 1 &mu;m)
 - EGM96 height (meters, accurate to 1 &mu;m)
 - EGM2008 height (meters, accurate to 1 &mu;m)
 .
The latitude and longitude are all multiples of 10<sup>&minus;6</sup>
deg and should be regarded as exact.  The geoid heights are computed
using the harmonic NGA synthesis programs, where the programs were
compiled (with gfortran) and run under Linux.  In the case of the
EGM2008 program, a <code>SAVE</code> statement needed to be added to
subroutine <code>latf</code>, in order for the program to compile
correctly with a stack-based compiler.  Similarly the EGM96 code
requires a <code>SAVE</code> statement in subroutine
<code>legfdn</code>.

<center>
Back to \ref other.  Forward to \ref gravity.  Up to \ref contents.
</center>
**********************************************************************/
/**
\page gravity Gravity models

<center>
Back to \ref geoid.  Forward to \ref magnetic.  Up to \ref contents.
</center>

%GeographicLib can compute the earth's gravitational field with an
earth gravity model using the GeographicLib::GravityModel and
GeographicLib::GravityCircle classes and with the
<a href="Gravity.1.html">Gravity</a> utility.  These models expand the
gravitational potential of the earth as sum of spherical harmonics.  The
models also specify a reference ellipsoid, relative to which geoid
heights and gravity disturbances are measured.

The supported models are
 - <b>egm84</b>, the
   <a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html">
   Earth Gravity Model 1984</a>, which includes terms up to degree 180.
 - <b>egm96</b>, the
   <a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html">
   Earth Gravity Model 1996</a>, which includes terms up to degree 360.
 - <b>egm2008</b>, the
   <a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008">
   Earth Gravity Model 2008</a>, which includes terms up to degree 2190.
 - <b>wgs84</b>, the
   <a href="http://earth-info.nga.mil/GandG/publications/tr8350.2/tr8350_2.html">
   WGS84 Reference Ellipsoid</a>.  This is just reproduces the normal
   gravitational field for the reference ellipsoid.  Usually
   GeographicLib::NormalGravity::WGS84 should be used instead.

See
 - W. A. Heiskanen and H. Moritz, Physical Geodesy (Freeman, San
   Francisco, 1967).
 .
for more information.

<b>Acknowledgment:</b> I would like to thank Mathieu Peyr&eacute;ga for
sharing EGM_Geoid_CalculatorClass from his Geo library with me.  His
implementation was the first I could easily understand and he and I
together worked through some of the issues with overflow and underflow
the occur while performing high-degree spherical harmonic sums.

Go to
 - \ref gravityinst
 - \ref gravityformat
 - \ref gravitynga
 - \ref gravitygeoid
 - \ref gravityatmos
 - \ref gravityparallel

\section gravityinst Installing the gravity models

These gravity models are available for download:
<center>
<table>
<caption>Available gravity models</caption>
<tr>
 <th rowspan="2">name         <th rowspan="2">max\n degree
 <th rowspan="2">size\n(kB)
 <th colspan="3"><center>Download Links (size, kB)</center></th>
<tr>
 <th>tar file</th>
 <th>Windows\n installer</th>
 <th>zip file</th>
<tr>
 <td>egm84
 <td><center>180</center>
 <td><center>27</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/gravity-distrib/egm84.tar.bz2">
 link</a> (26)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/gravity-distrib/egm84.exe">
 link</a> (55)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/gravity-distrib/egm84.zip">
 link</a> (26)</center>
<tr>
 <td>egm96
 <td><center>360</center>
 <td><center>2100</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/gravity-distrib/egm96.tar.bz2">
 link</a> (2100)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/gravity-distrib/egm96.exe">
 link</a> (2300)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/gravity-distrib/egm96.zip">
 link</a> (2100)</center>
<tr>
 <td>egm2008
 <td><center>2190</center>
 <td><center>76000</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/gravity-distrib/egm2008.tar.bz2">
 link</a> (75000)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/gravity-distrib/egm2008.exe">
 link</a> (72000)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/gravity-distrib/egm2008.zip">
 link</a> (73000)</center>
<tr>
 <td>wgs84
 <td><center>20</center>
 <td><center>1</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/gravity-distrib/wgs84.tar.bz2">
 link</a> (1)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/gravity-distrib/wgs84.exe">
 link</a> (30)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/gravity-distrib/wgs84.zip">
 link</a> (1)</center>
</table>
</center>
The "size" column is the size of the uncompressed data.

For Linux and Unix systems, %GeographicLib provides a shell script
geographiclib-get-gravity (typically installed in /usr/local/sbin)
which automates the process of downloading and installing the gravity
models.  For example
\verbatim
   geographiclib-get-gravity all  # to install egm84, egm96, egm2008, wgs84
   geographiclib-get-gravity -h   # for help
\endverbatim
This script should be run as a user with write access to the
installation directory, which is typically
/usr/local/share/GeographicLib (this can be overridden with the -p
flag), and the data will then be placed in the "gravity" subdirectory.

Windows users should download and run the Windows installers.  These
will prompt for an installation directory with the default being one of
\verbatim
   C:/Documents and Settings/All Users/Application Data/GeographicLib
   C:/ProgramData/GeographicLib
\endverbatim
(which you probably should not change) and the data is installed in the
"gravity" sub-directory.  (The second directory name is an alternate name
that Windows 7 for the "Application Data" directory.)

Otherwise download \e either the tar.bz2 file \e or the zip file (they
have the same contents).  To unpack these, run, for example
\verbatim
   mkdir -p /usr/local/share/GeographicLib
   tar xofjC egm96.tar.bz2 /usr/local/share/GeographicLib
   tar xofjC egm2008.tar.bz2 /usr/local/share/GeographicLib
   etc.
\endverbatim
and, again, the data will be placed in the "gravity" subdirectory.

However you install the gravity models, all the datasets should be
installed in the same directory.  GeographicLib::GravityModel and
<a href="Gravity.1.html">Gravity</a> uses a compile time default to
locate the datasets.  This is
- /usr/local/share/GeographicLib/gravity, for non-Windows systems
- C:/Documents and Settings/All Users/Application Data/GeographicLib/gravity,
  for Windows systems
.
consistent with the examples above.  This may be overridden at run-time
by defining the GRAVITY_PATH or the GEOGRAPHIC_DATA environment
variables; see GeographicLib::GravityModel::DefaultGravityPath() for
details.  Finally, the path may be set using the optional second
argument to the GeographicLib::GravityModel constructor or with the "-d"
flag to
<a href="Gravity.1.html">Gravity</a>.  Supplying the "-h" flag to
<a href="Gravity.1.html">Gravity</a> reports the default path for
gravity models for that utility.  The "-v" flag causes Gravity to report
the full path name of the data file it uses.

\section gravityformat The format of the gravity model files

The constructor for GeographicLib::GravityModel reads a file called
NAME.egm which specifies various properties for the gravity model.  It
then opens a binary file NAME.egm.cof to obtain the coefficients of the
spherical harmonic sum.

The first line of the .egm file must consist of "EGMF-v" where EGMF
stands for "Earth Gravity Model Format" and v is the version number of
the format (currently "1").

The rest of the File is read a line at a time.  A # character and
everything after it are discarded.  If the result is just white space it
is discarded.  The remaining lines are of the form "KEY WHITESPACE
VALUE".  In general, the KEY and the VALUE are case-sensitive.

GeographicLib::GravityModel only pays attention to the following
keywords
 - keywords that affect the field calculation, namely:
   - <b>ModelRadius</b> (required), the normalizing radius for the model
     in meters.
   - <b>ReferenceRadius</b> (required), the major radius \e a for the
     reference ellipsoid meters.
   - <b>ModelMass</b> (required), the mass constant \e GM for the model
     in meters<sup>3</sup>/seconds<sup>2</sup>.
   - <b>ReferenceMass</b> (required), the mass constant \e GM for the
     reference ellipsoid in meters<sup>3</sup>/seconds<sup>2</sup>.
   - <b>AngularVelocity</b> (required), the angular velocity &omega; for
     the model \e and the reference ellipsoid in rad
     seconds<sup>&minus;1</sup>.
   - <b>Flattening</b>, the flattening of the reference ellipsoid; this
     can be given a fraction, e.g., 1/298.257223563.  One of
     <b>Flattening</b> and <b>DynamicalFormFactor</b> is required.
   - <b>DynamicalFormFactor</b>, the dynamical form factor
     <i>J</i><sub>2</sub> for the reference ellipsoid.  One of
     <b>Flattening</b> and <b>DynamicalFormFactor</b> is required.
   - <b>HeightOffset</b> (default 0), the constant height offset
     (meters) added to obtain the geoid height.
   - <b>CorrectionMultiplier</b> (default 1), the multiplier for the
     "zeta-to-N" correction terms for the geoid height to convert them
     to meters.
   - <b>Normalization</b> (default full), the normalization used for the
     associated Legendre functions (full or schmidt).
   - <b>ID</b> (required), 8 printable characters which serve as a
     signature for the .egm.cof file (they must appear as the first 8
     bytes of this file).
   .
   The parameters <b>ModelRadius</b>, <b>ModelMass</b>, and
   <b>AngularVelocity</b> apply to the gravity model, while
   <b>ReferenceRadius</b>, <b>ReferenceMass</b>, <b>AngularVelocity</b>,
   and either <b>Flattening</b> or <b>DynamicalFormFactor</b>
   characterize the reference ellipsoid.  <b>AngularVelocity</b>
   (because it can be measured precisely) is the same for the gravity
   model and the reference ellipsoid.  <b>ModelRadius</b> is merely a
   scaling parameter for the gravity model and there's no requirement
   that it be close to the major radius of the earth (although that's
   typically how it is chosen).  <b>ModelMass</b> and
   <b>ReferenceMass</b> need not be the same and, indeed, they are
   slightly difference for egm2008.  As a result the disturbing
   potential \e T has a 1/\e r dependence at large distances.
 - keywords that store data that the user can query:
   - <b>Name</b>, the name of the model.
   - <b>Description</b>, a more descriptive name of the model.
   - <b>ReleaseDate</b>, when the model was created.
 - keywords that are examined to verify that their values are valid:
   - <b>ByteOrder</b> (default little), the order of bytes in the
     .egm.cof file.  Only little endian is supported at present.
 .
Other keywords are ignored.

The coefficient file NAME.egm.cof is a binary file in little endian
order.  The first 8 bytes of this file must match the ID given in
NAME.egm.  This is followed by 2 sets of spherical harmonic
coefficients.  The first of these gives the gravity potential and the
second gives the zeta-to-N corrections to the geoid height.  The format
for each set of coefficients is:
 - \e N, the maximum degree of the sum stored as a 4-byte signed integer.
   This must satisfy \e N &ge; &minus;1.
 - \e M, the maximum order of the sum stored as a 4-byte signed integer.
   This must satisfy \e N &ge; \e M &ge; &minus;1.
 - \e C<sub>\e nm</sub>, the coefficients of the cosine coefficients of
   the sum in column (i.e., \e m) major order.  There are (\e M + 1)
   (2\e N - \e M + 2) / 2 elements which are stored as IEEE doubles (8
   bytes).  For example for \e N = \e M = 3, there are 10 coefficients
   arranged as
   <i>C</i><sub>00</sub>,
   <i>C</i><sub>10</sub>,
   <i>C</i><sub>20</sub>,
   <i>C</i><sub>30</sub>,
   <i>C</i><sub>11</sub>,
   <i>C</i><sub>21</sub>,
   <i>C</i><sub>31</sub>,
   <i>C</i><sub>22</sub>,
   <i>C</i><sub>32</sub>,
   <i>C</i><sub>33</sub>.
 - \e S<sub>\e nm</sub>, the coefficients of the sine coefficients of
   the sum in column (i.e., \e m) major order starting at \e m = 1.
   There are \e M (2\e N - \e M + 1) / 2 elements which are stored as
   IEEE doubles (8 bytes).  For example for \e N = \e M = 3, there are 6
   coefficients arranged as
   <i>S</i><sub>11</sub>,
   <i>S</i><sub>21</sub>,
   <i>S</i><sub>31</sub>,
   <i>S</i><sub>22</sub>,
   <i>S</i><sub>32</sub>,
   <i>S</i><sub>33</sub>.
 .
Although the coefficient file is in little endian order, %GeographicLib
can read it on big endian machines.  It can only be read on machines
which store doubles in IEEE format.

As an illustration, here is egm2008.egm:
\verbatim
EGMF-1
# An Earth Gravity Model (Format 1) file.  For documentation on the
# format of this file see
# http://geographiclib.sf.net/html/gravity.html#gravityformat
Name            egm2008
Publisher       National Geospatial Intelligence Agency
Description     Earth Gravity Model 2008
URL             http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008
ReleaseDate     2008-06-01
ConversionDate  2011-11-19
DataVersion     1
ModelRadius     6378136.3
ModelMass       3986004.415e8
AngularVelocity 7292115e-11
ReferenceRadius 6378137
ReferenceMass   3986004.418e8
Flattening      1/298.257223563
HeightOffset    -0.41

# Gravitational and correction coefficients taken from
# EGM2008_to2190_TideFree and Zeta-to-N_to2160_egm2008 from
# the egm2008 distribution.
ID              EGM2008A
\endverbatim

\section gravitynga Comments on the NGA harmonic synthesis code

GeographicLib::GravityModel attempts to reproduce the results of NGA's
harmonic synthesis code for EGM2008, hsynth_WGS84.f.  Listed here are
issues that I encountered using the NGA code:
 -# A compiler which allocates local variables on the stack produces an
    executable with just returns NaNs.  The problem here is a missing
    <code>SAVE</code> statement in subroutine <code>latf</code>.
 -# Because the model and references masses for egm2008 differ (by about
    1 part in 10<sup>9</sup>), there should be a 1/\e r contribution to
    the disturbing potential \e T.  However, this term is set to zero in
    hsynth_WGS84 (effectively altering the normal potential).  This
    shifts the surface \e W = <i>U</i><sub>0</sub> outward by about 5 mm.
    Note too that the reference ellipsoid is no longer a level surface of
    the effective normal potential.
 -# Subroutine <code>radgrav</code> computes the ellipsoidal coordinate
    &beta; incorrectly.  This leads to small errors in the deflection
    of the vertical, &xi; and &eta;, when the height above the
    ellipsoid, \e h, is non-zero (about 10<sup>&minus;7</sup> arcsec at
    \e h = 400 km).
 -# There are several expressions which will return inaccurate results
    due to cancellation.  For example, subroutine <code>grs</code>
    computes the flattening using \e f = 1 - sqrt(1 -
    <i>e</i><sup>2</sup>).  Much better is to use \e f =
    <i>e</i><sup>2</sup>/(1 + sqrt(1 - <i>e</i><sup>2</sup>)).  The
    expressions for \e q and \e q' in <code>grs</code> and
    <code>radgrav</code> suffer from similar problems.  The resulting
    errors are tiny (about 50 pm in the geoid height); however, given
    that's there's essentially no cost to using more accurate
    expressions, it's preferable to do so.
 -# hsynth_WGS84 returns an "undefined" value for \e xi and \e eta at the
    poles.  Better would be to return the value obtained by taking the
    limit \e lat -> +/- 90&deg;.
 .
Issues 1--4 have been reported to the authors of hsynth_WGS84.
Issue 1 is peculiar to Fortran and is not encountered in C++ code and
GeographicLib::GravityModel corrects issues 3--5.  On issue 2,
GeographicLib::GravityModel neglects the 1/\e r term in \e T in
GeographicLib::GravityModel::GeoidHeight and
GeographicLib::GravityModel::SphericalAnomaly in order to produce
results which match NGA's for these quantities.  On the other hand,
GeographicLib::GravityModel::Disturbance and
GeographicLib::GravityModel::T <i>do</i> include this term.

\section gravitygeoid Details of the geoid height and anomaly calculations

Ideally, the geoid represents a surface of constant gravitational
potential which approximates mean sea level.  In reality some
approximations are taken in determining this surface.  The steps taking
by GeographicLib::GravityModel in computing the geoid height are
described here (in the context of EGM2008).  This mimics NGA's code
hsynth_WGS84 closely because most users of EGM2008 use the gridded data
generated by this code (e.g., GeographicLib::Geoid) and it is desirable
to use a consistent definition of the geoid height.
 - The model potential is band limited; the minimum wavelength that is
   represented is 360&deg;/2160 = 10' (i.e., about 10NM or
   18.5km).  The maximum degree for the spherical harmonic sum is 2190;
   however the model was derived using ellipsoidal harmonics of degree
   and order 2160 and the degree was increased to 2190 in order to
   capture the ellipsoidal harmonics faithfully with spherical
   harmonics.
 - The 1/\e r term is omitted from the \e T (this is issue 2 in \ref
   gravitynga).  This moves the equipotential surface by about 5mm.
 - The surface \e W = <i>U</i><sub>0</sub> is determined by Bruns'
   formula, which is roughly equivalent to a single iteration of
   Newton's method.  The RMS error in this approximation is about 1.5mm
   with a maximum error of about 10 mm.
 - The model potential is only valid above the earth's surface.  A
   correction therefore needs to be included where the geoid lies
   beneath the terrain.  This is NGA's "zeta-to-N" correction, which is
   represented by a spherical harmonic sum of degree and order 2160 (and
   so it misses short wavelength terrain variations).  In addition, it
   entails estimating the isostatic equilibrium of the earth's crust.
   The correction lies in the range [-5.05 m, 0.05 m], however for 99.9%
   of the earth's surface the correction is less than 10 mm in
   magnitude.
 - The resulting surface lies above the observed mean sea level,
   so -0.41m is added to the geoid height.  (Better would be to change
   the potential used to define the geoid; but this would only change
   the result by about 2mm.)
 .
A useful discussion of the problems with defining a geoid is given by
Dru A. Smith in
<a href="http://www.ngs.noaa.gov/PUBS_LIB/EGM96_GEOID_PAPER/egm96_geoid_paper.html">
There is no such thing as "The" EGM96 geoid: Subtle points on the use of
a global geopotential model</a>, IGeS Bulletin No. 8, International
Geoid Service, Milan, Italy, pp. 17--28 (1998).

GeographicLib::GravityModel::GeoidHeight reproduces the results of the
several NGA codes for harmonic synthesis with the following maximum
discrepancies:
 - egm84 = 1.1mm.  This is probably due to inconsistent parameters for the
   reference ellipsoid in the NGA code.  (In particular, the value of
   mass constant excludes the atmosphere; however, it's not clear
   whether the other parameters have been correspondingly adjusted.)
   Note that geoid heights predicted by egm84 differ from those of more
   recent gravity models by about 1 meter.
 - egm96 = 23nm.
 - egm2008 = 78pm.  After addressing some of the issues alluded to in
   issue 4 in \ref gravitynga, the maximum discrepancy becomes 23pm.

The formula for the gravity anomaly vector involves computing gravity
and normal gravity at two different points (with the displacement
between the points unknown <i>ab initio</i>).  Since the gravity anomaly
is already a small quantity it is sometimes acceptable to employ
approximations that change the quantities by \e O(\e f).  The NGA code
uses the spherical approximation described by Heiskanen and Moritz,
Sec. 2-14 and GeographicLib::GravityModel::SphericalAnomaly uses the
same approximation for compatibility.  In this approximation, the
gravity disturbance <b>delta</b> = <b>grad</b> \e T is calculated.
Here, \e T once again excludes the 1/\e r term (this is issue 2 in \ref
gravitynga and is consistent with the computation of the geoid height).
Note that <b>delta</b> compares the gravity and the normal gravity at
the \e same point; the gravity anomaly vector is then found by
estimating the gradient of the normal gravity in the limit that the
earth is spherically symmetric.  <b>delta</b> is expressed in \e
spherical coordinates as \e deltax, \e deltay, \e deltaz where, for
example, \e deltaz is the \e radial component of <b>delta</b> (not the
component perpendicular to the ellipsoid) and \e deltay is similarly
slightly different from the usual northerly component.  The components
of the anomaly are then given by
 - gravity anomaly, \e Dg01 = \e deltaz - 2<i>T</i>/\e R, where \e R
   distance to the center of the earth;
 - northerly component of the deflection of the vertical, \e xi = -
   <i>deltay</i>/\e gamma, where \e gamma is the magnitude of the normal
   gravity;
 - easterly component of the deflection of the vertical, \e eta = -
   <i>deltax</i>/\e gamma.
 .
GeographicLib::NormalGravity computes the normal gravity accurately and
avoids issue 3 of \ref gravitynga.  Thus while
GeographicLib::GravityModel::SphericalAnomaly reproduces the results for
\e xi and \e eta at \e h = 0, there is a slight discrepancy if \e h is
non-zero.

\section gravityatmos The effect of the mass of the atmosphere

All of the supported models use WGS84 for the reference ellipsoid.  This
has (see
<a href="http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf">
TR8350.2</a>, table 3.1)
 - \e a = 6378137 m
 - \e f = 1/298.257223563
 - &omega; = 7292115 &times; 10<sup>&minus;11</sup> rad
   s<sup>&minus;1</sup>
 - \e GM = 3986004.418 &times; 10<sup>8</sup> m<sup>3</sup>
   s<sup>&minus;2</sup>.
 .
The value of \e GM includes the mass of the atmosphere and so strictly
only applies above the earth's atmosphere.  Near the surface of the
earth, the value of \e g will be less (in absolute value) than the value
predicted by these models by about &delta;\e g = (4&pi;<i>G</i>/\e
g) \e A = 8.552 &times; 10<sup>&minus;11</sup> \e A m<sup>2</sup>/kg,
where \e G is the gravitational constant, \e g is the earth's gravity,
and \e A is the pressure of the atmosphere.  At sea level we have \e A =
101.3 kPa, and &delta;\e g = 8.7 &times; 10<sup>&minus;6</sup> m
s<sup>&minus;2</sup>, approximately.  (In other words the effect is
about 1&nbsp;part in a million; by way of comparison, buoyancy effects
are about 100 times larger.)

\section gravityparallel Geoid heights on a multi-processor system

The egm2008 model includes many terms (over 2 million spherical
harmonics).  For that reason computations using this model may be slow;
for example it takes about 78 ms to compute the geoid height at a single
point.  There are two ways to speed up this computation:
 - Use a GeographicLib::GravityCircle to compute the geoid height at
   several points on a circle of latitude.  This reduces the cost per
   point to about 92 &mu;s (a reduction by a factor of over 800).
 - Compute the values on several circles of latitude in parallel.  One
   of the simplest ways of doing this is with
   <a href="http://openmp.org"> OpenMP</a>; on an 8-processor system,
   this can speed up the computation by another factor of 8.
 .
Both of these techniques are illustrated by the following code,
which computes a table of geoid heights on
a regular grid and writes on the result in a
<a href="http://vdatum.noaa.gov/dev/gtx_info.html#dev_gtx_binary">.gtx</a>
file.  On an 8-processor Intel 2.66 GHz machine using OpenMP
(-DHAVE_OPENMP=1), it takes about 40 minutes of elapsed time to compute
the geoid height for EGM2008 on a 1' gride.  (Without these
optimizations, the computation would have taken about 200 days!)
\include GeoidToGTX.cpp

This example, <code>examples/GeoidToGTX.cpp</code>, is built if cmake
is configured with <code>-D GEOGRAPHICLIB_EXAMPLES=ON</code>.  cmake
will add in support for OpenMP, if it is available.

<center>
Back to \ref geoid.  Forward to \ref magnetic.  Up to \ref contents.
</center>
**********************************************************************/
/**
\page magnetic Magnetic models

<center>
Back to \ref gravity.  Forward to \ref geodesic.  Up to \ref contents.
</center>

%GeographicLib can compute the earth's magnetic field by a magnetic
model using the GeographicLib::MagneticModel and
GeographicLib::MagneticCircle classes and with the
<a href="MagneticField.1.html">MagneticField</a> utility.  These models
expand the internal magnetic potential of the earth as sum of spherical
harmonics.  They neglect magnetic fields due to the ionosphere, the
magnetosphere, nearby magnetized materials, electric machinery, etc.
Users of GeographicLib::MagneticModel are advised to read the
<a href="http://www.ngdc.noaa.gov/IAGA/vmod/igrfhw.html">"Health
Warning"</a> this is provided with igrf11.  Although the advice is
specific to igrf11, many of the comments apply to all magnetic field
models.

The supported models are
 - <b>wmm2010</b>, the
   <a href="http://ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml"> World
   Magnetic Model 2010</a>, which approximates the main magnetic field
   for the period 2010--2015.
 - <b>igrf11</b>, the
   <a href="http://ngdc.noaa.gov/IAGA/vmod/igrf.html">International
   Geomagnetic Reference Field (11th generation)</a>, which approximates
   the main magnetic field for the period 1900--2015.
 - <b>emm2010</b>, the
   <a href="http://ngdc.noaa.gov/geomag/EMM/index.html"> Enhanced
   Magnetic Model 2010</a>, which approximates the main and crustal
   magnetic fields for the period 2010--2015.

Go to
 - \ref magneticinst
 - \ref magneticformat

\section magneticinst Installing the magnetic field models

These magnetic models are available for download:
<center>
<table>
<caption>Available magnetic models</caption>
<tr>
 <th rowspan="2">name         <th rowspan="2">max\n degree
 <th rowspan="2">time\n interval
 <th rowspan="2">size\n(kB)
 <th colspan="3"><center>Download Links (size, kB)</center></th>
<tr>
 <th>tar file</th>
 <th>Windows\n installer</th>
 <th>zip file</th>
<tr>
 <td>wmm2010
 <td><center>12</center>
 <td><center>2010--2015</center>
 <td><center>3</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/magnetic-distrib/wmm2010.tar.bz2">
 link</a> (2)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/magnetic-distrib/wmm2010.exe">
 link</a> (300)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/magnetic-distrib/wmm2010.zip">
 link</a> (2)</center>
<tr>
 <td>igrf11
 <td><center>13</center>
 <td><center>1900--2015</center>
 <td><center>25</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/magnetic-distrib/igrf11.tar.bz2">
 link</a> (7)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/magnetic-distrib/igrf11.exe">
 link</a> (310)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/magnetic-distrib/igrf11.zip">
 link</a> (8)</center>
<tr>
 <td>emm2010
 <td><center>739</center>
 <td><center>2010--2015</center>
 <td><center>4400</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/magnetic-distrib/emm2010.tar.bz2">
 link</a> (3700)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/magnetic-distrib/emm2010.exe">
 link</a> (3000)</center>
 <td><center>
 <a href="http://sf.net/projects/geographiclib/files/magnetic-distrib/emm2010.zip">
 link</a> (4100)</center>
</table>
</center>
The "size" column is the size of the uncompressed data.

For Linux and Unix systems, %GeographicLib provides a shell script
geographiclib-get-magnetic (typically installed in /usr/local/sbin)
which automates the process of downloading and installing the magnetic
models.  For example
\verbatim
   geographiclib-get-magnetic all  # to install wmm2010, igrf11, emm2010
   geographiclib-get-magnetic -h   # for help
\endverbatim
This script should be run as a user with write access to the
installation directory, which is typically
/usr/local/share/GeographicLib (this can be overridden with the -p
flag), and the data will then be placed in the "magnetic" subdirectory.

Windows users should download and run the Windows installers.  These
will prompt for an installation directory with the default being one of
\verbatim
   C:/Documents and Settings/All Users/Application Data/GeographicLib
   C:/ProgramData/GeographicLib
\endverbatim
(which you probably should not change) and the data is installed in the
"magnetic" sub-directory.  (The second directory name is an alternate name
that Windows 7 for the "Application Data" directory.)

Otherwise download \e either the tar.bz2 file \e or the zip file (they
have the same contents).  To unpack these, run, for example
\verbatim
   mkdir -p /usr/local/share/GeographicLib
   tar xofjC wmm2010.tar.bz2 /usr/local/share/GeographicLib
   tar xofjC emm2010.tar.bz2 /usr/local/share/GeographicLib
   etc.
\endverbatim
and, again, the data will be placed in the "magnetic" subdirectory.

However you install the magnetic models, all the datasets should be
installed in the same directory.  GeographicLib::MagneticModel and
<a href="MagneticField.1.html">MagneticField</a> uses a compile time
default to locate the datasets.  This is
- /usr/local/share/GeographicLib/magnetic, for non-Windows systems
- C:/Documents and Settings/All Users/Application Data/GeographicLib/magnetic,
  for Windows systems
.
consistent with the examples above.  This may be overridden at run-time
by defining the MAGNETIC_PATH  or the GEOGRAPHIC_DATA environment
variables; see GeographicLib::MagneticModel::DefaultMagneticPath() for
details.  Finally, the path may be set using the optional second
argument to the GeographicLib::MagneticModel constructor or with the
"-d" flag to <a href="MagneticField.1.html">MagneticField</a>.
Supplying the "-h" flag to
<a href="MagneticField.1.html">MagneticField</a> reports the default
path for magnetic models for that utility.  The "-v" flag causes
MagneticField to report the full path name of the data file it uses.

\section magneticformat The format of the magnetic model files

The constructor for GeographicLib::MagneticModel reads a file called
NAME.wmm which specifies various properties for the magnetic model.  It
then opens a binary file NAME.wmm.cof to obtain the coefficients of the
spherical harmonic sum.

The first line of the .wmm file must consist of "WMMF-v" where WMMF
stands for "World Magnetic Model Format" and v is the version number of
the format (currently "1").

The rest of the File is read a line at a time.  A # character and
everything after it are discarded.  If the result is just white space it
is discarded.  The remaining lines are of the form "KEY WHITESPACE
VALUE".  In general, the KEY and the VALUE are case-sensitive.

GeographicLib::MagneticModel only pays attention to the following
keywords
 - keywords that affect the field calculation, namely:
   - <b>Radius</b> (required), the normalizing radius of the model in
     meters.
   - <b>NumModels</b> (default 1), the number of models.  WMM2010
     consists of a single model giving the magnetic field and its time
     variation at 2010.  IGRF11 consists of 23 models for 1900 thru 2010
     at 5 year intervals.  The time variation is given only for the last
     model to allow extrapolation beyond 2010.  For dates prior to 2010,
     linear interpolation is used.
   - <b>Epoch</b> (required), the time origin (in fractional years) for
     the first model.
   - <b>DeltaEpoch</b> (default 1), the interval between models in years
     (only relevant for NumModels > 1).
   - <b>Normalization</b> (default schmidt), the normalization used for
     the associated Legendre functions (schmidt or full).
   - <b>ID</b> (required), 8 printable characters which serve as a
     signature for the .wmm.cof file (they must appear as the first 8
     bytes of this file).
 - keywords that store data that the user can query:
   - <b>Name</b>, the name of the model.
   - <b>Description</b>, a more descriptive name of the model.
   - <b>ReleaseDate</b>, when the model was created.
   - <b>MinTime</b>, the minimum date at which the model should be used.
   - <b>MaxTime</b>, the maximum date at which the model should be used.
   - <b>MinHeight</b>, the minimum height above the ellipsoid for which
     the model should be used.
   - <b>MaxHeight</b>, the maximum height above the ellipsoid for which
     the model should be used.
   .
   GeographicLib::MagneticModel does not enforce the restrictions
   implied by last four quantities.  However,
   <a href="MagneticField.1.html">MagneticField</a> issues a warning if
   these limits are exceeded.
 - keywords that are examined to verify that their values are valid:
   - <b>Type</b> (default linear), the type of the model.  "linear"
     means that the time variation is piece-wise linear (either using
     interpolation between the field at two dates or using the field and
     its first derivative with respect to time).  This is the only type
     of model supported at present.
   - <b>ByteOrder</b> (default little), the order of bytes in the
     .wmm.cof file.  Only little endian is supported at present.
 .
Other keywords are ignored.

The coefficient file NAME.wmm.cof is a binary file in little endian
order.  The first 8 bytes of this file must match the ID given in
NAME.wmm.  This is followed by NumModels + 1 sets of spherical harmonic
coefficients.  The first NumModels of these model the magnetic field at
Epoch + \e i * DeltaEpoch for 0 &le; \e i < NumModels.  The last set of
coefficients model the rate of change of the magnetic field at Epoch +
(NumModels &minus; 1) * DeltaEpoch.  The format for each set of coefficients
is:
 - \e N, the maximum degree of the sum stored as a 4-byte signed integer.
   This must satisfy \e N &ge; &minus;1.
 - \e M, the maximum order of the sum stored as a 4-byte signed integer.
   This must satisfy \e N &ge; \e M &ge; &minus;1.
 - \e C<sub>\e nm</sub>, the coefficients of the cosine coefficients of
   the sum in column (i.e., \e m) major order.  There are (\e M + 1)
   (2\e N &minus; \e M + 2) / 2 elements which are stored as IEEE doubles (8
   bytes).  For example for \e N = \e M = 3, there are 10 coefficients
   arranged as
   <i>C</i><sub>00</sub>,
   <i>C</i><sub>10</sub>,
   <i>C</i><sub>20</sub>,
   <i>C</i><sub>30</sub>,
   <i>C</i><sub>11</sub>,
   <i>C</i><sub>21</sub>,
   <i>C</i><sub>31</sub>,
   <i>C</i><sub>22</sub>,
   <i>C</i><sub>32</sub>,
   <i>C</i><sub>33</sub>.
 - \e S<sub>\e nm</sub>, the coefficients of the sine coefficients of
   the sum in column (i.e., \e m) major order starting at \e m = 1.
   There are \e M (2\e N &minus; \e M + 1) / 2 elements which are stored as
   IEEE doubles (8 bytes).  For example for \e N = \e M = 3, there are 6
   coefficients arranged as
   <i>S</i><sub>11</sub>,
   <i>S</i><sub>21</sub>,
   <i>S</i><sub>31</sub>,
   <i>S</i><sub>22</sub>,
   <i>S</i><sub>32</sub>,
   <i>S</i><sub>33</sub>.
 .
Although the coefficient file is in little endian order, %GeographicLib
can read it on big endian machines.  It can only be read on machines
which store doubles in IEEE format.

As an illustration, here is igrf11.wmm:
\verbatim
WMMF-1
# A World Magnetic Model (Format 1) file.  For documentation on the
# format of this file see
# http://geographiclib.sf.net/html/magnetic.html#magneticformat
Name            igrf11
Description     International Geomagnetic Reference Field 11th Generation
URL             http://ngdc.noaa.gov/IAGA/vmod/igrf.html
Publisher       National Oceanic and Atmospheric Administration
ReleaseDate     2009-12-15
DataCutOff      2009-10-01
ConversionDate  2011-11-04
DataVersion     1
Radius          6371200
NumModels       23
Epoch           1900
DeltaEpoch      5
MinTime         1900
MaxTime         2015
MinHeight       -1000
MaxHeight       600000

# The coefficients are stored in a file obtained by appending ".cof" to
# the name of this file.  The coefficients were obtained from IGRF11.COF
# in the geomag70 distribution.
ID              IGRF11-A
\endverbatim

<center>
Back to \ref gravity.  Forward to \ref geodesic.  Up to \ref contents.
</center>
**********************************************************************/
/**
\page geodesic Geodesics on an ellipsoid of revolution

<center>
Back to \ref magnetic.  Forward to \ref triaxial.  Up to \ref contents.
</center>

GeographicLib::Geodesic and GeographicLib::GeodesicLine provide accurate
solutions to the direct and inverse geodesic problems.  The
<a href="GeodSolve.1.html">GeodSolve</a> utility provides an interface
to these classes.  GeographicLib::AzimuthalEquidistant implements the
azimuthal equidistant projection in terms of geodesics.
GeographicLib::CassiniSoldner implements a transverse cylindrical
equidistant projection in terms of geodesics.  The <a
href="GeodesicProj.1.html">GeodesicProj</a> utility provides an
interface to these projections.

The algorithms used by GeographicLib::Geodesic and
GeographicLib::GeodesicLine are based on a Taylor expansion of the
geodesic integrals valid when the flattening \e f is small.
GeographicLib::GeodesicExact and GeographicLib::GeodesicLineExact
evaluate the integrals exactly (in terms of incomplete elliptic
integrals).  For the WGS84 ellipsoid, the series solutions are about
2--3 times faster and 2--3 times more accurate (because it's easier to
control round-off errors with series solutions); thus
GeographicLib::Geodesic and GeographicLib::GeodesicLine are
recommended for most geodetic applications.  However, in applications
where the absolute value of \e f is greater than about 0.02, the exact
classes should be used.

Go to
 - \ref testgeod
 - \ref geodseries
 - \ref geodellip
 - \ref meridian
 .
For some background information on geodesics on triaxial ellipsoids, see
\ref triaxial.

References
 - F. W. Bessel,
   <a href="http://dx.doi.org/10.1002/asna.201011352">The calculation
   of longitude and latitude from geodesic measurements (1825)</a>,
   Astron. Nachr. 331(8), 852-861 (2010);
   translated by C. F. F. Karney and R. E. Deakin.  Preprint:
   <a href="http://arxiv.org/abs/0908.1824">arXiv:0908.1824</a>.
 - F. R. Helmert,
   <a href="http://geographiclib.sf.net/geodesic-papers/helmert80-en.pdf">
   Mathematical and Physical Theories of Higher Geodesy, Part 1 (1880)</a>,
   Aeronautical Chart and Information Center (St. Louis, 1964),
   Chaps. 5--7.
 - J. Danielsen,
   The Area under the Geodesic,
   Survey Review 30(232), 61--66 (1989).
 - C. F. F. Karney,
   <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
   Algorithms for geodesics</a>,
   J. Geodesy 87(1), 43--55 (Jan. 2013);
   DOI: <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
   10.1007/s00190-012-0578-z</a>;
   addenda: <a href="http://geographiclib.sf.net/geod-addenda.html">
   geod-addenda.html</a>;
   resource page:
   <a href="http://geographiclib.sf.net/geod.html"> geod.html</a>.
 - A collection of some papers on geodesics is available at
   http://geographiclib.sourceforge.net/geodesic-papers/biblio.html

\section testgeod Test data for geodesics

A test set a geodesics is available at
 - <a href="http://sf.net/projects/geographiclib/files/testdata/GeodTest.dat.gz">
   GeodTest.dat.gz</a>
 .
This is about 39 MB (compressed).  This consists of a set of geodesics
for the WGS84 ellipsoid.  A subset of this (consisting of 1/50 of the
members &mdash; about 690 kB, compressed) is available at
 - <a href="http://sf.net/projects/geographiclib/files/testdata/GeodTest-short.dat.gz">
   GeodTest-short.dat.gz</a>

Each line of the test set gives 10 space delimited numbers
 - latitude for point 1, \e lat1 (degrees, exact)
 - longitude for point 1, \e lon1 (degrees, always 0)
 - azimuth for point 1, \e azi1 (clockwise from north in degrees, exact)
 - latitude for point 2, \e lat2 (degrees, accurate to
   10<sup>&minus;18</sup> deg)
 - longitude for point 2, \e lon2 (degrees, accurate to
   10<sup>&minus;18</sup> deg)
 - azimuth for point 2, \e azi2 (degrees, accurate to
   10<sup>&minus;18</sup> deg)
 - geodesic distance from point 1 to point 2, \e s12 (meters, exact)
 - arc distance on the auxiliary sphere, \e a12 (degrees, accurate to
   10<sup>&minus;18</sup> deg)
 - reduced length of the geodesic, \e m12 (meters, accurate to 0.1 pm)
 - the area under the geodesic, \e S12 (m<sup>2</sup>, accurate to
   1 mm<sup>2</sup>)
 .
These are computed using as direct geodesic calculations with the given
\e lat1, \e lon1, \e azi1, and \e s12.  The distance \e s12 always
corresponds to an arc length \e a12 &le; 180&deg;, so the given
geodesics give the shortest paths from point 1 to point 2.  For
simplicity and without loss of generality, \e lat1 is chosen in
[0&deg;, 90&deg;], \e lon1 is taken to be zero, \e azi1 is
chosen in [0&deg;, 180&deg;].  Furthermore, \e lat1 and \e
azi1 are taken to be multiples of 10<sup>&minus;12</sup> deg and \e s12
is a multiple of 0.1 &mu;m in [0 m, 20003931.4586254 m].  This results \e
lon2 in [0&deg;, 180&deg;] and \e azi2 in [0&deg;, 180&deg;].

The direct calculation uses an expansion of the geodesic equations
accurate to <i>f</i><sup>30</sup> (approximately 1 part in 10<sup>50</sup>)
and is computed with with
<a href="http://en.wikipedia.org/wiki/Maxima_(software)">Maxima</a>'s
bfloats and fpprec set to 100 (so the errors in the data are probably
1/2 of the values quoted above).

The contents of the file are as follows:
 - 100000 entries randomly distributed
 - 50000 entries which are nearly antipodal
 - 50000 entries with short distances
 - 50000 entries with one end near a pole
 - 50000 entries with both ends near opposite poles
 - 50000 entries which are nearly meridional
 - 50000 entries which are nearly equatorial
 - 50000 entries running between vertices (\e azi1 = \e azi2 = 90&deg;)
 - 50000 entries ending close to vertices
 .
(a total of 500000 entries).  The values for \e s12 for the geodesics
running between vertices are truncated to a multiple of 0.1 pm and this
is used to determine point 2.

This data can be fed to the <a href="GeodSolve.1.html">GeodSolve</a>
utility as follows
 - Direct from point 1:
\verbatim
  gunzip -c GeodTest.dat.gz | cut -d' ' -f1,2,3,7 | ./GeodSolve
\endverbatim
   This should yield columns 4, 5, 6, and 9 of the test set.
 - Direct from point 2:
\verbatim
  gunzip -c GeodTest.dat.gz | cut -d' ' -f4,5,6,7 |
  sed "s/ \([^ ]*$\)/ -\1/" | ./GeodSolve
\endverbatim
   (The sed command negates the distance.)  This should yield columns 1,
   2, and 3, and the negative of column 9 of the test set.
 - Inverse between points 1 and 2:
\verbatim
  gunzip -c GeodTest.dat.gz | cut -d' ' -f1,2,4,5 | ./GeodSolve -i
\endverbatim
   This should yield columns 3, 6, 7, and 9 of the test set.
 .
Add, e.g., "-p 6", to the call to GeodSolve to change the precision of
the output.  Adding "-f" causes GeodSolve to print 12 fields specifying
the geodesic; these include the 10 fields in the test set plus the
geodesic scales \e M12 and \e M21 which are inserted between \e m12 and
\e S12.

Code for computing arbitrarily accurate geodesics in maxima is available
in <a href="geodesic.mac"> geodesic.mac</a> (this depends on
<a href="ellint.mac"> ellint.mac</a> and uses the series computed by
<a href="geod.mac"> geod.mac</a>).  This solve both the direct and
inverse geodesic problems and offers the ability to solve the problems
either using series expansions (similar to GeographicLib::Geodesic) or
in terms of elliptic integrals (similar to
GeographicLib::GeodesicExact).

\section geodseries Expansions for geodesics

We give here the series expansions for the various geodesic integrals
valid to order <i>f</i><sup>10</sup>.  In this release of the code, we
use a 6th-order expansions.  This is sufficient to maintain accuracy
for doubles for the SRMmax ellipsoid (\e a = 6400 km, \e f = 1/150).
However, the preprocessor macro GEOGRAPHICLIB_GEODESIC_ORDER can be
used to select any order up to 8.  (If using long doubles, with a
64-bit fraction, the default order is 7.)  The series expanded to
order <i>f</i><sup>30</sup> are given in <a href="geodseries30.html">
geodseries30.html</a>.

In the formulas below ^ indicates exponentiation (\e f^3 =
<i>f</i><sup>3</sup>) and / indicates real division (3/5 = 0.6).  The
equations need to be converted to Horner form, but are here left in
expanded form so that they can be easily truncated to lower order.
These expansions were obtained using the Maxima code,
<a href="geod.mac">geod.mac</a>.

In the expansions below, we have
 - \f$ \alpha \f$ is the azimuth
 - \f$ \alpha_0 \f$ is the azimuth at the equator crossing
 - \f$ \lambda \f$ is the longitude measured from the equator crossing
 - \f$ \sigma \f$ is the spherical arc length
 - \f$ \omega = \tan^{-1}(\sin\alpha_0\tan\sigma) \f$ is the spherical longitude
 - \f$ a \f$ is the equatorial radius
 - \f$ b \f$ is the polar semi-axis
 - \f$ f \f$ is the flattening
 - \f$ e^2 = f(2 - f) \f$
 - \f$ e'^2 = e^2/(1-e^2) \f$
 - \f$ k^2 = e'^2 \cos^2\alpha_0 = 4 \epsilon / (1 - \epsilon)^2 \f$
 - \f$ n = f / (2 - f) \f$
 - \f$ c^2 = a^2/2 + b^2/2 (\tanh^{-1}e)/e \f$
 - \e ep2 = \f$ e'^2 \f$
 - \e k2 = \f$ k^2 \f$
 - \e eps = \f$ \epsilon \f$

The formula for distance is
\f[
 \frac sb = I_1(\sigma)
\f]
where
\f[
\begin{aligned}
  I_1(\sigma) &= A_1\bigl(\sigma + B_1(\sigma)\bigr) \\
  B_1(\sigma) &= \sum_{j=1} C_{1j} \sin 2j\sigma
\end{aligned}
\f]
and
\verbatim
A1 = (1 + 1/4 * eps^2
        + 1/64 * eps^4
        + 1/256 * eps^6
        + 25/16384 * eps^8
        + 49/65536 * eps^10) / (1 - eps);
\endverbatim
\verbatim
C1[1] = - 1/2 * eps
        + 3/16 * eps^3
        - 1/32 * eps^5
        + 19/2048 * eps^7
        - 3/4096 * eps^9;
C1[2] = - 1/16 * eps^2
        + 1/32 * eps^4
        - 9/2048 * eps^6
        + 7/4096 * eps^8
        + 1/65536 * eps^10;
C1[3] = - 1/48 * eps^3
        + 3/256 * eps^5
        - 3/2048 * eps^7
        + 17/24576 * eps^9;
C1[4] = - 5/512 * eps^4
        + 3/512 * eps^6
        - 11/16384 * eps^8
        + 3/8192 * eps^10;
C1[5] = - 7/1280 * eps^5
        + 7/2048 * eps^7
        - 3/8192 * eps^9;
C1[6] = - 7/2048 * eps^6
        + 9/4096 * eps^8
        - 117/524288 * eps^10;
C1[7] = - 33/14336 * eps^7
        + 99/65536 * eps^9;
C1[8] = - 429/262144 * eps^8
        + 143/131072 * eps^10;
C1[9] = - 715/589824 * eps^9;
C1[10] = - 2431/2621440 * eps^10;
\endverbatim

The function \f$ \tau(\sigma) = s/(b A_1) = \sigma + B_1(\sigma) \f$
may be inverted by series reversion giving
\f[
  \sigma(\tau) = \tau + \sum_{j=1} C'_{1j} \sin 2j\sigma
\f]
where
\verbatim
C1'[1] = + 1/2 * eps
         - 9/32 * eps^3
         + 205/1536 * eps^5
         - 4879/73728 * eps^7
         + 9039/327680 * eps^9;
C1'[2] = + 5/16 * eps^2
         - 37/96 * eps^4
         + 1335/4096 * eps^6
         - 86171/368640 * eps^8
         + 4119073/28311552 * eps^10;
C1'[3] = + 29/96 * eps^3
         - 75/128 * eps^5
         + 2901/4096 * eps^7
         - 443327/655360 * eps^9;
C1'[4] = + 539/1536 * eps^4
         - 2391/2560 * eps^6
         + 1082857/737280 * eps^8
         - 2722891/1548288 * eps^10;
C1'[5] = + 3467/7680 * eps^5
         - 28223/18432 * eps^7
         + 1361343/458752 * eps^9;
C1'[6] = + 38081/61440 * eps^6
         - 733437/286720 * eps^8
         + 10820079/1835008 * eps^10;
C1'[7] = + 459485/516096 * eps^7
         - 709743/163840 * eps^9;
C1'[8] = + 109167851/82575360 * eps^8
         - 550835669/74317824 * eps^10;
C1'[9] = + 83141299/41287680 * eps^9;
C1'[10] = + 9303339907/2972712960 * eps^10;
\endverbatim

The reduced length is given by
\f[
\begin{aligned}
  \frac mb &= \sqrt{1 + k^2 \sin^2\sigma_2} \cos\sigma_1 \sin\sigma_2 \\
  &\quad {}-\sqrt{1 + k^2 \sin^2\sigma_1} \sin\sigma_1 \cos\sigma_2 \\
  &\quad {}-\cos\sigma_1 \cos\sigma_2 \bigl(J(\sigma_2) - J(\sigma_1)\bigr)
\end{aligned}
\f]
where
\f[
\begin{aligned}
 J(\sigma) &= I_1(\sigma) - I_2(\sigma) \\
 I_2(\sigma) &= A_2\bigl(\sigma + B_2(\sigma)\bigr) \\
 B_2(\sigma) &= \sum_{j=1} C_{2j} \sin 2j\sigma
\end{aligned}
\f]
\verbatim
A2 = (1 + 1/4 * eps^2
        + 9/64 * eps^4
        + 25/256 * eps^6
        + 1225/16384 * eps^8
        + 3969/65536 * eps^10) * (1 - eps);
\endverbatim
\verbatim
C2[1] = + 1/2 * eps
        + 1/16 * eps^3
        + 1/32 * eps^5
        + 41/2048 * eps^7
        + 59/4096 * eps^9;
C2[2] = + 3/16 * eps^2
        + 1/32 * eps^4
        + 35/2048 * eps^6
        + 47/4096 * eps^8
        + 557/65536 * eps^10;
C2[3] = + 5/48 * eps^3
        + 5/256 * eps^5
        + 23/2048 * eps^7
        + 191/24576 * eps^9;
C2[4] = + 35/512 * eps^4
        + 7/512 * eps^6
        + 133/16384 * eps^8
        + 47/8192 * eps^10;
C2[5] = + 63/1280 * eps^5
        + 21/2048 * eps^7
        + 51/8192 * eps^9;
C2[6] = + 77/2048 * eps^6
        + 33/4096 * eps^8
        + 2607/524288 * eps^10;
C2[7] = + 429/14336 * eps^7
        + 429/65536 * eps^9;
C2[8] = + 6435/262144 * eps^8
        + 715/131072 * eps^10;
C2[9] = + 12155/589824 * eps^9;
C2[10] = + 46189/2621440 * eps^10;
\endverbatim

The longitude is given in terms of the spherical longitude by
\f[
\lambda = \omega - f \sin\alpha_0 I_3(\sigma)
\f]
where
\f[
\begin{aligned}
 I_3(\sigma) &= A_3\bigl(\sigma + B_3(\sigma)\bigr) \\
 B_3(\sigma) &= \sum_{j=1} C_{3j} \sin 2j\sigma
\end{aligned}
\f]
and
\verbatim
A3 = 1 - (1/2 - 1/2*n) * eps
       - (1/4 + 1/8*n - 3/8*n^2) * eps^2
       - (1/16 + 3/16*n + 1/16*n^2 - 5/16*n^3) * eps^3
       - (3/64 + 1/32*n + 5/32*n^2 + 5/128*n^3 - 35/128*n^4) * eps^4
       - (3/128 + 5/128*n + 5/256*n^2 + 35/256*n^3 + 7/256*n^4) * eps^5
       - (5/256 + 15/1024*n + 35/1024*n^2 + 7/512*n^3) * eps^6
       - (25/2048 + 35/2048*n + 21/2048*n^2) * eps^7
       - (175/16384 + 35/4096*n) * eps^8
       - 245/32768 * eps^9;
\endverbatim
\verbatim
C3[1] = + (1/4 - 1/4*n) * eps
        + (1/8 - 1/8*n^2) * eps^2
        + (3/64 + 3/64*n - 1/64*n^2 - 5/64*n^3) * eps^3
        + (5/128 + 1/64*n + 1/64*n^2 - 1/64*n^3 - 7/128*n^4) * eps^4
        + (3/128 + 11/512*n + 3/512*n^2 + 1/256*n^3 - 7/512*n^4) * eps^5
        + (21/1024 + 5/512*n + 13/1024*n^2 + 1/512*n^3) * eps^6
        + (243/16384 + 189/16384*n + 83/16384*n^2) * eps^7
        + (435/32768 + 109/16384*n) * eps^8
        + 345/32768 * eps^9;
C3[2] = + (1/16 - 3/32*n + 1/32*n^2) * eps^2
        + (3/64 - 1/32*n - 3/64*n^2 + 1/32*n^3) * eps^3
        + (3/128 + 1/128*n - 9/256*n^2 - 3/128*n^3 + 7/256*n^4) * eps^4
        + (5/256 + 1/256*n - 1/128*n^2 - 7/256*n^3 - 3/256*n^4) * eps^5
        + (27/2048 + 69/8192*n - 39/8192*n^2 - 47/4096*n^3) * eps^6
        + (187/16384 + 39/8192*n + 31/16384*n^2) * eps^7
        + (287/32768 + 47/8192*n) * eps^8
        + 255/32768 * eps^9;
C3[3] = + (5/192 - 3/64*n + 5/192*n^2 - 1/192*n^3) * eps^3
        + (3/128 - 5/192*n - 1/64*n^2 + 5/192*n^3 - 1/128*n^4) * eps^4
        + (7/512 - 1/384*n - 77/3072*n^2 + 5/3072*n^3 + 65/3072*n^4) * eps^5
        + (3/256 - 1/1024*n - 71/6144*n^2 - 47/3072*n^3) * eps^6
        + (139/16384 + 143/49152*n - 383/49152*n^2) * eps^7
        + (243/32768 + 95/49152*n) * eps^8
        + 581/98304 * eps^9;
C3[4] = + (7/512 - 7/256*n + 5/256*n^2 - 7/1024*n^3 + 1/1024*n^4) * eps^4
        + (7/512 - 5/256*n - 7/2048*n^2 + 9/512*n^3 - 21/2048*n^4) * eps^5
        + (9/1024 - 43/8192*n - 129/8192*n^2 + 39/4096*n^3) * eps^6
        + (127/16384 - 23/8192*n - 165/16384*n^2) * eps^7
        + (193/32768 + 3/8192*n) * eps^8
        + 171/32768 * eps^9;
C3[5] = + (21/2560 - 9/512*n + 15/1024*n^2 - 7/1024*n^3 + 9/5120*n^4) * eps^5
        + (9/1024 - 15/1024*n + 3/2048*n^2 + 57/5120*n^3) * eps^6
        + (99/16384 - 91/16384*n - 781/81920*n^2) * eps^7
        + (179/32768 - 55/16384*n) * eps^8
        + 141/32768 * eps^9;
C3[6] = + (11/2048 - 99/8192*n + 275/24576*n^2 - 77/12288*n^3) * eps^6
        + (99/16384 - 275/24576*n + 55/16384*n^2) * eps^7
        + (143/32768 - 253/49152*n) * eps^8
        + 33/8192 * eps^9;
C3[7] = + (429/114688 - 143/16384*n + 143/16384*n^2) * eps^7
        + (143/32768 - 143/16384*n) * eps^8
        + 429/131072 * eps^9;
C3[8] = + (715/262144 - 429/65536*n) * eps^8
        + 429/131072 * eps^9;
C3[9] = + 2431/1179648 * eps^9;
\endverbatim

The formula for area between the geodesic and the equator is given in
Sec. 6 of <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
Algorithms for geodesics</a> in terms of \e S,
\f[
S = c^2 \alpha + e^2 a^2 \cos\alpha_0 \sin\alpha_0 I_4(\sigma)
\f]
where
\f[
I_4(\sigma) = \sum_{j=0} C_{4j} \cos(2j+1)\sigma
\f]
In the paper, this was expanded in \f$ e'^2 \f$ and \f$ k^2 \f$.
However, the series converges faster for eccentric ellipsoids if the
expansion is in \f$ n \f$ and \f$ \epsilon \f$.  The series to order
\f$ f^{10} \f$ becomes
\verbatim
C4[0] = + (2/3 - 4/15*n + 8/105*n^2 + 4/315*n^3 + 16/3465*n^4 + 20/9009*n^5 + 8/6435*n^6 + 28/36465*n^7 + 32/62985*n^8 + 4/11305*n^9)
        - (1/5 - 16/35*n + 32/105*n^2 - 16/385*n^3 - 64/15015*n^4 - 16/15015*n^5 - 32/85085*n^6 - 112/692835*n^7 - 128/1616615*n^8) * eps
        - (2/105 + 32/315*n - 1088/3465*n^2 + 1184/5005*n^3 - 128/3465*n^4 - 3232/765765*n^5 - 1856/1616615*n^6 - 6304/14549535*n^7) * eps^2
        + (11/315 - 368/3465*n - 32/6435*n^2 + 976/4095*n^3 - 154048/765765*n^4 + 368/11115*n^5 + 5216/1322685*n^6) * eps^3
        + (4/1155 + 1088/45045*n - 128/1287*n^2 + 64/3927*n^3 + 2877184/14549535*n^4 - 370112/2078505*n^5) * eps^4
        + (97/15015 - 464/45045*n + 4192/153153*n^2 - 88240/969969*n^3 + 31168/1322685*n^4) * eps^5
        + (10/9009 + 4192/765765*n - 188096/14549535*n^2 + 23392/855855*n^3) * eps^6
        + (193/85085 - 6832/2078505*n + 106976/14549535*n^2) * eps^7
        + (632/1322685 + 3456/1616615*n) * eps^8
        + 107/101745 * eps^9;
C4[1] = + (1/45 - 16/315*n + 32/945*n^2 - 16/3465*n^3 - 64/135135*n^4 - 16/135135*n^5 - 32/765765*n^6 - 112/6235515*n^7 - 128/14549535*n^8) * eps
        - (2/105 - 64/945*n + 128/1485*n^2 - 1984/45045*n^3 + 256/45045*n^4 + 64/109395*n^5 + 128/855855*n^6 + 2368/43648605*n^7) * eps^2
        - (1/105 - 16/2079*n - 5792/135135*n^2 + 3568/45045*n^3 - 103744/2297295*n^4 + 264464/43648605*n^5 + 544/855855*n^6) * eps^3
        + (4/1155 - 2944/135135*n + 256/9009*n^2 + 17536/765765*n^3 - 3053056/43648605*n^4 + 1923968/43648605*n^5) * eps^4
        + (1/9009 + 16/19305*n - 2656/153153*n^2 + 65072/2078505*n^3 + 526912/43648605*n^4) * eps^5
        + (10/9009 - 1472/459459*n + 106112/43648605*n^2 - 204352/14549535*n^3) * eps^6
        + (349/2297295 + 28144/43648605*n - 32288/8729721*n^2) * eps^7
        + (632/1322685 - 44288/43648605*n) * eps^8
        + 43/479655 * eps^9;
C4[2] = + (4/525 - 32/1575*n + 64/3465*n^2 - 32/5005*n^3 + 128/225225*n^4 + 32/765765*n^5 + 64/8083075*n^6 + 32/14549535*n^7) * eps^2
        - (8/1575 - 128/5775*n + 256/6825*n^2 - 6784/225225*n^3 + 4608/425425*n^4 - 128/124355*n^5 - 5888/72747675*n^6) * eps^3
        - (8/1925 - 1856/225225*n - 128/17325*n^2 + 42176/1276275*n^3 - 2434816/72747675*n^4 + 195136/14549535*n^5) * eps^4
        + (8/10725 - 128/17325*n + 64256/3828825*n^2 - 128/25935*n^3 - 266752/10392525*n^4) * eps^5
        - (4/25025 + 928/3828825*n + 292288/72747675*n^2 - 106528/6613425*n^3) * eps^6
        + (464/1276275 - 17152/10392525*n + 83456/72747675*n^2) * eps^7
        + (1168/72747675 + 128/1865325*n) * eps^8
        + 208/1119195 * eps^9;
C4[3] = + (8/2205 - 256/24255*n + 512/45045*n^2 - 256/45045*n^3 + 1024/765765*n^4 - 256/2909907*n^5 - 512/101846745*n^6) * eps^3
        - (16/8085 - 1024/105105*n + 2048/105105*n^2 - 1024/51051*n^3 + 4096/373065*n^4 - 1024/357357*n^5) * eps^4
        - (136/63063 - 256/45045*n + 512/1072071*n^2 + 494336/33948915*n^3 - 44032/1996995*n^4) * eps^5
        + (64/315315 - 16384/5360355*n + 966656/101846745*n^2 - 868352/101846745*n^3) * eps^6
        - (16/97461 + 14848/101846745*n + 74752/101846745*n^2) * eps^7
        + (5024/33948915 - 96256/101846745*n) * eps^8
        - 1744/101846745 * eps^9;
C4[4] = + (64/31185 - 512/81081*n + 1024/135135*n^2 - 512/109395*n^3 + 2048/1247103*n^4 - 2560/8729721*n^5) * eps^4
        - (128/135135 - 2048/405405*n + 77824/6891885*n^2 - 198656/14549535*n^3 + 8192/855855*n^4) * eps^5
        - (512/405405 - 2048/530145*n + 299008/130945815*n^2 + 280576/43648605*n^3) * eps^6
        + (128/2297295 - 2048/1438965*n + 241664/43648605*n^2) * eps^7
        - (17536/130945815 + 1024/43648605*n) * eps^8
        + 2944/43648605 * eps^9;
C4[5] = + (128/99099 - 2048/495495*n + 4096/765765*n^2 - 6144/1616615*n^3 + 8192/4849845*n^4) * eps^5
        - (256/495495 - 8192/2807805*n + 376832/53348295*n^2 - 8192/855855*n^3) * eps^6
        - (6784/8423415 - 432128/160044885*n + 397312/160044885*n^2) * eps^7
        + (512/53348295 - 16384/22863555*n) * eps^8
        - 16768/160044885 * eps^9;
C4[6] = + (512/585585 - 4096/1422135*n + 8192/2078505*n^2 - 4096/1322685*n^3) * eps^6
        - (1024/3318315 - 16384/9006855*n + 98304/21015995*n^2) * eps^7
        - (103424/189143955 - 8192/4203199*n) * eps^8
        - 1024/189143955 * eps^9;
C4[7] = + (1024/1640925 - 65536/31177575*n + 131072/43648605*n^2) * eps^7
        - (2048/10392525 - 262144/218243025*n) * eps^8
        - 84992/218243025 * eps^9;
C4[8] = + (16384/35334585 - 131072/82447365*n) * eps^8
        - 32768/247342095 * eps^9;
C4[9] = + 32768/92147055 * eps^9;
\endverbatim

\section geodellip Geodesics in terms of elliptic integrals

GeographicLib::GeodesicExact and GeographicLib::GeodesicLineExact solve
the geodesic problem using elliptic integrals.  The formulation of
geodesic in terms of incomplete elliptic integrals is given in
 - C. F. F. Karney,
   <a href="http://arxiv.org/abs/1102.1215v1">Geodesics
   on an ellipsoid of revolution</a>,
   Feb. 2011; preprint
   <a href="http://arxiv.org/abs/1102.1215v1">arxiv:1102.1215v1</a>.
 .
It is most convenient to use the form derived for a prolate ellipsoid in
Appendix D.  For an oblate ellipsoid this results in elliptic integrals
with an imaginary modulus.  However, the integrals themselves are real
and the algorithms used to compute the elliptic integrals handles the
case of an imaginary modulus using real arithmetic.

The key relations used by %GeographicLib are
\f[
 \begin{aligned}
  \frac sb &= E(\sigma, ik), \\
  \lambda &= (1 - f) \sin\alpha_0  G(\sigma, \cos^2\alpha_0, ik) \\
          &= \chi
           - \frac{e'^2}{\sqrt{1+e'^2}}\sin\alpha_0 H(\sigma, -e'^2, ik), \\
  J(\sigma) &= k^2 D(\sigma, ik),
 \end{aligned}
\f]
where \f$ \chi \f$ is a modified spherical longitude given by
\f[
\tan\chi = \sqrt{\frac{1+e'^2}{1+k^2\sin^2\sigma}}\tan\omega,
\f]
and
\f[
 \begin{aligned}
 D(\phi,k) &= \int_0^\phi
 \frac{\sin^2\theta}{\sqrt{1 - k^2\sin^2\theta}}\,d\theta\\
 &=\frac{F(\phi, k) - E(\phi, k)}{k^2},\\
 G(\phi,\alpha^2,k) &= \int_0^\phi
 \frac{\sqrt{1 - k^2\sin^2\theta}}{1 - \alpha^2\sin^2\theta}\,d\theta\\
 &=\frac{k^2}{\alpha^2}F(\phi, k)
 +\biggl(1-\frac{k^2}{\alpha^2}\biggr)\Pi(\phi, \alpha^2, k),\\
 H(\phi, \alpha^2, k)
 &= \int_0^\phi
   \frac{\cos^2\theta}{(1-\alpha^2\sin^2\theta)\sqrt{1-k^2\sin^2\theta}}
   \,d\theta \\
 &=
 \frac1{\alpha^2} F(\phi, k) +
      \biggl(1 - \frac1{\alpha^2}\biggr) \Pi(\phi, \alpha^2, k),
 \end{aligned}
\f]
and \f$F(\phi, k)\f$, \f$E(\phi, k)\f$, \f$D(\phi, k)\f$, and
\f$\Pi(\phi, \alpha^2, k)\f$, are incomplete elliptic integrals (see
http://dlmf.nist.gov/19.2.ii).  The formula for \f$ s \f$ and the
first expression for \f$ \lambda \f$ are given by Legendre (1811) and
are the most common representation of geodesics in terms of elliptic
integrals.  The second (equivalent) expression for \f$ \lambda \f$,
which was given by Cayley (1870), is useful in that the elliptic
integral is relegated to a small correction term.  This form allows
the longitude to be computed more accurately and is used in
%GeographicLib.  (The equivalence of the two expressions for \f$
\lambda \f$ follows from http://dlmf.nist.gov/19.7.E8.)

Nominally, GeographicLib::GeodesicExact and
GeographicLib::GeodesicLineExact will give "exact" results for any value
of the flattening.  However, the geographic latitude is a distorted
measure of distance from the equator with very eccentric ellipsoids and
this introducing an irreducible representational error in the algorithms
in this case.  It is therefore recommended to restrict the use of these
classes to \e b/\e a &isin; [0.01, 100] or \e f &isin; [-99, 0.99].
Note that GeographicLib::GeodesicExact still uses a series expansion for
the area \e S12.  However the series is taken out to 30th order and
gives accurate results for \e b/\e a &isin; [1/2, 2]; the accuracy is
about 8 decimal digits for \e b/\e a &isin; [1/4, 4].  Additional work
planned for this aspect of the geodesic problem:
- formulate the area integral \e S12 in terms of elliptic integrals;
- generate accurate test geodesics for highly eccentric ellipsoids so
  that the roundoff errors can be quantified.

Thomas (1952) and Rollins (2010) use a different independent variable
for geodesics, \f$\theta\f$ instead of \f$\sigma\f$, where \f$
\tan\theta = \sqrt{1 + k^2} \tan\sigma \f$.  The corresponding
expressions for \f$ s \f$ and \f$ \lambda \f$ are given here for
completeness:
\f[
\begin{aligned}
\frac sb &= \sqrt{1-k'^2} \Pi(\theta, k'^2, k'), \\
\lambda &= (1-f) \sqrt{1-k'^2} \sin\alpha_0 \Pi(\theta, k'^2/e^2, k'),
\end{aligned}
\f]
where \f$ k' = k/\sqrt{1 + k^2} \f$.  The expression for \f$ s \f$
can be written in terms of elliptic integrals of the second kind and
Cayley's technique can be used to subtract out the leading order
behavior of \f$ \lambda \f$ to give
\f[
\begin{aligned}
\frac sb &=\frac1{\sqrt{1-k'^2}}
  \biggl( E(\theta, k') -
    \frac{k'^2 \sin\theta \cos\theta}{\sqrt{1-k'^2\sin^2\theta}} \biggr), \\
\lambda &= \psi + (1-f) \sqrt{1-k'^2} \sin\alpha_0
\bigl( F(\theta, k') - \Pi(\theta, e^2, k') \bigr),
\end{aligned}
\f]
where
\f[
\begin{aligned}
\tan\psi &= \sqrt{\frac{1+k^2\sin^2\sigma}{1+e'^2}}\tan\omega \\
         &= \sqrt{\frac{1-e^2}{1+k^2\cos^2\theta}}\sin\alpha_0\tan\theta.
\end{aligned}
\f]
The tangents of the three "longitude-like" angles are in geometric
progression, \f$ \tan\chi/\tan\omega = \tan\omega/\tan\psi \f$.

\section meridian Parameters for the meridian

The formulas for \f$ s \f$ given in the previous section are the same as
those for the distance along a meridian for an ellipsoid with equatorial
radius \f$ a \sqrt{1 - e^2 \sin^2\alpha_0} \f$ and polar semi-axis \f$ b
\f$.  Here is a list of possible ways of expressing the meridian
distance in terms of elliptic integrals using the notation:
- \f$ a \f$, equatorial axis,
- \f$ b \f$, polar axis,
- \f$ e = \sqrt{(a^2 - b^2)/a^2} \f$, eccentricity,
- \f$ e' = \sqrt{(a^2 - b^2)/b^2} \f$, second eccentricity,
- \f$ \phi = \mathrm{am}(u, e) \f$, the geographic latitude,
- \f$ \phi' = \mathrm{am}(v', ie') = \pi/2 - \phi \f$,
  the geographic colatitude,
- \f$ \beta = \mathrm{am}(v, ie') \f$, the parametric latitude
  (\f$ \tan^2\beta = (1 - e^2) \tan^2\phi \f$),
- \f$ \beta' = \mathrm{am}(u', e) = \pi/2 - \beta \f$,
  the parametric colatitude,
- \f$ M \f$, the length of a quarter meridian (equator to pole),
- \f$ y \f$, the distance along the meridian (measured from the equator).
- \f$ y' = M -y \f$, the distance along the meridian (measured from the pole).
.
The eccentricities \f$ (e, e') \f$ are real (resp. imaginary) for
oblate (resp. prolate) ellipsoids.  The elliptic variables \f$(u,
u')\f$ and \f$(v, v')\f$ are defined by
- \f$ u = F(\phi, e) ,\quad u' = F(\beta', e) \f$
- \f$ v = F(\beta, ie') ,\quad v' = F(\phi', ie') \f$,
.
and are linearly related by
- \f$ u + u' = K(e) ,\quad v + v' = K(ie') \f$
- \f$ v = \sqrt{1-e^2} u ,\quad u = \sqrt{1+e'^2} v \f$.
.
The cartesian coordinates for the meridian \f$ (x, y) \f$ are given by
\f[
\begin{aligned}
  x &= a \cos\beta = a \cos\phi / \sqrt{1 - e^2 \sin^2\phi} \\
    &= a \sin\beta' = (a^2/b) \sin\phi' / \sqrt{1 + e'^2 \sin^2\phi'} \\
    &= a \,\mathrm{cn}(v, ie) = a \,\mathrm{cd}(u, e) \\
    &= a \,\mathrm{sn}(u', e) = (a^2/b) \,\mathrm{sd}(v', ie'),
\end{aligned}
\f]
\f[
\begin{aligned}
  z &= b \sin\beta = (b^2/a) \sin\phi / \sqrt{1 - e^2 \sin^2\phi} \\
    &= b \cos\beta' = b \cos\phi' / \sqrt{1 + e'^2 \sin^2\phi'} \\
    &= b \,\mathrm{sn}(v, ie) = (b^2/a) \,\mathrm{sd}(u, e)  \\
    &= b \,\mathrm{cn}(u', e) = b \,\mathrm{cd}(v', ie').
\end{aligned}
\f]
The distance along the meridian can be expressed variously as
\f[
\begin{aligned}
  y  &=   b   \int \sqrt{1 + e'^2  \sin^2\beta}\, d\beta
     =   b E(\beta, ie') \\
     &= \frac{b^2}a \int \frac1{(1 - e^2 \sin^2\phi)^{3/2}}\, d\phi
      = \frac{b^2}a \Pi(\phi, e^2, e) \\
     &=   a \biggl(E(\phi, e) -
        \frac{e^2\sin\phi \cos\phi}{\sqrt{1 - e^2\sin^2\phi}}\biggr) \\
     &=  b \int \mathrm{dn}^2(v, ie')\, dv
      = \frac{b^2}a \int \mathrm{nd}^2(u, e)\, du
      = \cal E(v, ie'),
\end{aligned}
\f]
\f[
\begin{aligned}
  y' &=   a   \int \sqrt{1 - e^2  \sin^2\beta'}\, d\beta'
     =   a E(\beta', e) \\
     &= \frac{a^2}b \int \frac1{(1 + e'^2 \sin^2\phi')^{3/2}}\, d\phi'
      = \frac{a^2}b \Pi(\phi', -e'^2, ie') \\
     &=  b \biggl(E(\phi', ie') +
        \frac{e'^2\sin\phi' \cos\phi'}{\sqrt{1 + e'^2\sin^2\phi'}}\biggr) \\
     &=  a \int \mathrm{dn}^2(u', e)\, du'
      = \frac{a^2}b \int \mathrm{nd}^2(v', ie')\, dv'
      = \cal E(u', e),
\end{aligned}
\f]
with the quarter meridian distance given by
\f[
  M  = aE(e) = bE(ie') = (b^2/a)\Pi(e^2,e) = (a^2/b)\Pi(-e'^2,ie').
\f]
(Here \f$ E, F, \Pi \f$ are elliptic integrals defined in
http://dlmf.nist.gov/19.2.ii.  \f$ \cal E, \mathrm{am}, \mathrm{sn},
\mathrm{cn}, \mathrm{sd}, \mathrm{cd}, \mathrm{dn}, \mathrm{nd} \f$ are
Jacobi elliptic functions defined in http://dlmf.nist.gov/22.2 and
http://dlmf.nist.gov/22.16.)

There are several considerations in the choice of independent variable
for evaluate the meridian distance
- The use of an imaginary modulus (namely, \f$ ie' \f$, above) is of no
  practical concern.  The integrals are real in this case and modern
  methods (%GeographicLib uses the method given in
  http://dlmf.nist.gov/19.36.i) for computing integrals handles this
  case using just real arithmetic.
- If the "natural" origin is the equator, choose one of \f$ \phi, \beta,
  u, v \f$ (this might be preferred in geodesy).  If it's the pole,
  choose one of the complementary quantities \f$ \phi', \beta', u', v'
  \f$ (this might be preferred by mathematicians).
- Applying these formulas to the geodesic problems, \f$ \beta \f$
  becomes the arc length, \f$ \sigma \f$, on the auxiliary sphere.  This
  is the traditional method of solution used by Legendre (1806), Oriani
  (1806), Bessel (1825), Helmert (1880), Rainsford (1955), Thomas
  (1970), Vincenty (1975), Rapp (1993), and so on.  Many of the
  solutions in terms of elliptic functions use one of the elliptic
  variables (\f$ u \f$ or \f$ v \f$), see, for example, Jacobi (1855),
  Halphen (1888), Forsyth (1896).  In the context of geodesics \f$
  \phi \f$ becomes Thomas' variable \f$ \theta \f$; this is used by
  Thomas (1952) and Rollins (2010) in their formulation of the
  geodesic problem (see the previous section).
- For highly eccentric ellipsoids the variation of the meridian with
  respect to \f$ \beta \f$ is considerably "better behaved" than other
  choices (see the figure below).  The choice of \f$ \phi \f$ is
  probably a poor one in this case.
.
%GeographicLib uses the geodesic generalization of
\f$ y = b E(\beta, ie') \f$, namely \f$ s = b E(\sigma, ik) \f$.  See
\ref geodellip.

<img src="meridian-measures.png" width=576 height=596
     alt="comparison of meridian measures">

<center>
Back to \ref magnetic.  Forward to \ref triaxial.  Up to \ref contents.
</center>
**********************************************************************/
/**
\page triaxial Geodesics on a triaxial ellipsoid

<center>
Back to \ref geodesic.  Forward to \ref transversemercator.  Up to
\ref contents.
</center>

Jacobi (1839) showed that the problem of geodesics on a triaxial
ellipsoids (with 3 unequal axes) can be reduced by quadrature.  Despite
this, the detailed behavior of the geodesics is not very well known.  In
this section, I briefly give Jacobi's solution and illustrate the
behavior of the geodesics and outline an algorithm for the solution of
the inverse problem.  The figures in this section are relatively small
in order to allow them to be displayed quickly.  Vector versions of
these figures are available in
<a href="http://geographiclib.sourceforge.net/triaxial-figs.pdf"
type="application/pdf"> triaxial-figs.pdf</a>.  This allows you to
magnify the figures to show the details more clearly.

Go to
 - \ref triaxial-coords
 - \ref triaxial-jacobi
 - \ref triaxial-survey
 - \ref triaxial-stab
 - \ref triaxial-inverse

<b>NOTES</b>
 -# A triaxial ellipsoid approximates the earth only slightly better
    than an ellipsoid of revolution.  If you are really considering
    measuring distances on the earth using a triaxial ellipsoid, you
    should also be worrying about the shape of the geoid, which
    essentially makes the geodesic problem a hopeless mess; see, for
    example, <a href="http://arxiv.org/abs/1112.3231"> Waters
    (2011)</a>.
 -# There is nothing new in this section.  It is just an exercise in
    exploring Jacobi's solution.  My interest here is in generating long
    geodesics with the correct long-time behavior.  Arnold gives a
    nice qualitative description of the solution in <i>Mathematical
    Methods of Classical Mechanics</i> (2nd edition, Springer, 1989),
    pp. 264--266.
 -# Possible reasons this problem might, nevertheless, be of interest
    are:
    - It is the first example of a dynamical system which has a
      non-trivial constant of motion.  As such, Jacobi's paper generated
      a lot of excitement and was followed by many papers elaborating
      his solution.  In particular, the unstable behavior of one of the
      closed geodesics of the ellipsoid, is an early example of a system
      with a positive Lyapunov exponent (one of the essential
      ingredients for chaotic behavior in dynamical systems).
    - Knowledge of ellipsoidal coordinates (used by Jacobi) might be
      useful in other areas of geodesy.
    - Geodesics which pass through the pole on an ellipsoid of revolution
      represent a degenerate class (they are all closed and all pass
      through the opposite pole).  It is of interest to see how this
      degeneracy is broken with a surface with a more general shape.
 -# My interest in this problem was piqued by Jean-Marc Baillard.  I put
    him onto Jacobi's solution without having looked at it in detail
    myself; and he quickly implemented the solution for an HP-41
    calculator(!) which is posted
    <a href="http://hp41programs.yolasite.com/geod3axial.php"> here</a>.
 -# I do not give full citations of the papers here.  You can find these
    in the
    <a href="http://geographiclib.sourceforge.net/geodesic-papers/biblio.html">
    Online Geodesic Bibliography</a>; this includes links to online
    versions of the papers.
 -# An alternative to exploring geodesics using Jacobi's solution is to
    integrate the equations for the geodesics directly.  This is the
    approach taken by
    <a href="http://www.math.harvard.edu/~knill/caustic/">
    Oliver Knill and Michael Teodorescu</a>.  However it is difficult to
    ensure that the long time behavior is correctly modeled with such an
    approach.
 -# At this point, I have no plans to add the solution of triaxial
    geodesic problem to %GeographicLib.
 -# If you only want to learn about geodesics on a biaxial ellipsoid (an
    ellipsoid of revolution), then see \ref geodesic or the paper
    - C. F. F. Karney,
      <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
      Algorithms for geodesics</a>,
      J. Geodesy 87(1), 43--55 (Jan. 2013);
      DOI: <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
      10.1007/s00190-012-0578-z</a>;
      addenda: <a href="http://geographiclib.sf.net/geod-addenda.html">
      geod-addenda.html</a>.

\section triaxial-coords Triaxial coordinate systems

Consider the ellipsoid defined by
\f[
  f = \frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1,
\f]
where, without loss of generality, \f$ a \ge b \ge c \gt 0\f$.  A
point on the surface is specified by a latitude and longitude.  The \e
geographical latitude and longitude \f$(\phi, \lambda)\f$ are defined by
\f[
 \frac{\nabla f}{\left| \nabla f\right|} = \left(
\begin{array}{c} \cos\phi \cos\lambda \\ \cos\phi \sin\lambda \\ \sin\phi
\end{array}\right).
\f]
The \e parametric latitude and longitude \f$(\phi', \lambda')\f$ are
defined by
\f[
\begin{aligned}
 x &= a \cos\phi' \cos\lambda', \\
 y &= b \cos\phi' \sin\lambda', \\
 z &= c \sin\phi'. \\
\end{aligned}
\f]
Jacobi employed the \e ellipsoidal latitude and longitude \f$(\beta,
\omega)\f$ defined by
\f[
\begin{aligned}
  x &= a \cos\omega
      \frac{\sqrt{a^2 - b^2\sin^2\beta - c^2\cos^2\beta}}
           {\sqrt{a^2 - c^2}}, \\
  y &= b \cos\beta \sin\omega, \\
  z &= c \sin\beta
      \frac{\sqrt{a^2\sin^2\omega + b^2\cos^2\omega - c^2}}
           {\sqrt{a^2 - c^2}}.
\end{aligned}
\f]
Grid lines of constant \f$\beta\f$ and \f$\omega\f$ are given in Fig. 1.

<center>
<img src="triaxial-grid.png"
width=419 height=397 alt="Geodesic grid on a triaxial ellipsoid">
Fig. 1
</center>\n
Fig. 1:
The ellipsoidal grid.  The blue (resp. green) lines are lines of constant
\f$\beta\f$ (resp. \f$\omega\f$); the grid spacing is 10&deg;.  Also
shown in red are two of the principal sections of the ellipsoid, defined
by \f$x = 0\f$ and \f$z = 0\f$.  The third principal section, \f$y =
0\f$, is covered by the lines \f$\beta = \pm 90^\circ\f$ and \f$\omega =
90^\circ \pm 90^\circ\f$.  These lines meet at four umbilical points (two
of which are visible in this figure) where the principal radii of
curvature are equal.  The parameters of the ellipsoid are \f$a =
1.01\f$, \f$b = 1\f$, \f$c = 0.8\f$, and it is viewed in an orthographic
projection from a point above \f$\phi = 40^\circ\f$, \f$\lambda =
30^\circ\f$.  These parameters were chosen to accentuate the ellipsoidal
effects on geodesics (relative to those on the earth) while still
allowing the connection to an ellipsoid of revolution to be made.

The grid lines of the ellipsoid coordinates are "lines of curvature" on
the ellipsoid, i.e., they are parallel to the direction of principal
curvature (Monge, 1796).  They are also intersections of the ellipsoid
with confocal systems of hyperboloids of one and two sheets (Dupin,
1813).  Finally they are geodesic ellipses and hyperbolas defined using
two adjacent umbilical points.  For example, the lines of constant
\f$\beta\f$ in Fig. 1 can be generated with the familiar string
construction for ellipses with the ends of the string pinned to the two
umbilical points.

Inter-conversions between these three latitudes and longitudes and the
cartesian coordinates are simple algebraic exercises.

\section triaxial-jacobi Jacobi's solution

Solving the geodesic problem for an ellipsoid of revolution is, from the
mathematical point of view, trivial; because of symmetry, geodesics have
a constant of the motion (analogous to the angular momentum) which was
found by Clairaut (1733).  By 1806 (with the work of Legendre, Oriani,
et al.), there was a complete understanding of the qualitative behavior
of geodesics on an ellipsoid of revolution.

On the other hand, geodesics on a triaxial ellipsoid have no obvious
constant of the motion and thus represented a challenging "unsolved"
problem in the first half of the nineteenth century.  Jacobi discovered
that the geodesic equations are separable if they are expressed in
ellipsoidal coordinates.  You can get an idea of the importance Jacobi
attached to his discovery from the
<a href="http://books.google.com/books?id=_09tAAAAMAAJ&pg=PA385">
letter</a> he wrote to his friend and neighbor Bessel:
<blockquote> The day before yesterday, I reduced to quadrature the
problem of geodesic lines on an <i>ellipsoid with three unequal
axes</i>.  They are the simplest formulas in the world, Abelian
integrals, which become the well known elliptic integrals if 2 axes are
set equal.\n
K&ouml;nigsberg, 28th Dec. '38.
</blockquote>

On the same day he wrote a similar letter to the editor of Compte Rendus
and his result was published in J. Crelle in (1839) with a French
translation (from German) appearing in J. Liouville in (1841).

Here is the solution, exactly as given by Jacobi
<a href="http://books.google.com/books?id=Rh8GAAAAYAAJ&pg=PA268"> here</a>
(with minor changes in notation):
\f[
\begin{aligned}
\delta &= \int \frac
{\sqrt{b^2\sin^2\beta + c^2\cos^2\beta}\,d\beta}
{\sqrt{a^2 - b^2\sin^2\beta - c^2\cos^2\beta}
 \sqrt{(b^2-c^2)\cos^2\beta - \gamma}}\\
&\quad -
\int \frac
{\sqrt{a^2\sin^2\omega + b^2\cos^2\omega}\,d\omega}
{\sqrt{a^2\sin^2\omega + b^2\cos^2\omega - c^2}
 \sqrt{(a^2-b^2)\sin^2\omega + \gamma}}
\end{aligned}
\f]
As Jacobi notes &quot;a function of the angle \f$\beta\f$ equals a
function of the angle \f$\omega\f$.  These two functions are just
Abelian integrals...&quot; Two constants \f$\delta\f$ and \f$\gamma\f$
appear in the solution.  Typically \f$\delta\f$ is zero if the lower
limits of the integrals are taken to be the starting point of the geodesic
and the direction of the geodesics is determined by \f$\gamma\f$.
However for geodesics that start at an umbilical points, we have \f$\gamma
= 0\f$ and \f$\delta\f$ determines the direction at the umbilical point.
Incidentally the constant \f$\gamma\f$ may be expressed as
\f[
\gamma = (b^2-c^2)\cos^2\beta\sin^2\alpha-(a^2-b^2)\sin^2\omega\cos^2\alpha
\f]
where \f$\alpha\f$ is the angle the geodesic makes with lines of
constant \f$\omega\f$. In the limit \f$a\rightarrow b\f$, this reduces
to \f$\cos\beta\sin\alpha = \text{const.}\f$, the familiar Clairaut
relation.  A nice derivation of Jacobi's result is given by Darboux
(1894) <a href="http://books.google.com/books?id=hGMSAAAAIAAJ&pg=PA9">
&sect;&sect;583--584</a> where he gives the solution found by Liouville
(1846) for general quadratic surfaces.  In this formulation, the
distance along the geodesic, \f$s\f$, is also found using
\f[
\begin{aligned}
\frac{ds}{(b^2-c^2)\cos^2\beta + (a^2-b^2)\sin^2\omega}
&= \frac
{\sqrt{b^2\sin^2\beta + c^2\cos^2\beta}\,d\beta}
{\sqrt{a^2 - b^2\sin^2\beta - c^2\cos^2\beta}
 \sqrt{(b^2-c^2)\cos^2\beta - \gamma}}\\
&= \frac
{\sqrt{a^2\sin^2\omega + b^2\cos^2\omega}\,d\omega}
{\sqrt{a^2\sin^2\omega + b^2\cos^2\omega - c^2}
 \sqrt{(a^2-b^2)\sin^2\omega + \gamma}}
\end{aligned}
\f]
An alternative expression for the distance is
\f[
\begin{aligned}
ds
&= \frac
{\sqrt{b^2\sin^2\beta + c^2\cos^2\beta}
 \sqrt{(b^2-c^2)\cos^2\beta - \gamma}\,d\beta}
{\sqrt{a^2 - b^2\sin^2\beta - c^2\cos^2\beta}}\\
&\quad {}+ \frac
{\sqrt{a^2\sin^2\omega + b^2\cos^2\omega}
 \sqrt{(a^2-b^2)\sin^2\omega + \gamma}\,d\omega}
{\sqrt{a^2\sin^2\omega + b^2\cos^2\omega - c^2}}
\end{aligned}
\f]

Jacobi's solution is a convenient way to compute geodesics on an
ellipsoid.  Care must be taken with the signs of the square roots (which
are determined by the initial azimuth of the geodesic).  Also if
\f$\gamma \gt 0\f$ (resp. \f$\gamma \lt 0\f$), then the \f$\beta\f$
(resp. \f$\omega\f$) integrand diverges.  The integrand can be
transformed into a finite one by a change of variable, e.g.,
\f$\sin\beta = \sin\sigma \sqrt{1 - \gamma/(b^2-c^2)}\f$.  The resulting
integrals are periodic, so the behavior of an arbitrarily long geodesic
is entirely captured by tabulating the integrals over a single period.

The situation is more complicated if \f$\gamma = 0\f$ (corresponding to
umbilical geodesics).  Both integrands have simple poles at the umbilical
points.  However, this behavior may be subtracted from the integrands to
yield (for example) the sum of a term involving
\f$\tanh^{-1}\sin\beta\f$ and a finite integral.  Since both integrals
contain similar logarithmic singularities they can be equated (thus
fixing the ratio \f$\cos\beta/\sin\omega\f$ at the umbilical point) and
connection formulas can be found which allow the geodesic to be followed
through the umbilical point.  The study of umbilical geodesics was of
special interest to a group of Irish mathematicians in the 1840's and
1850's, including Michael and William Roberts (twins!), Hart, Graves,
and Salmon.

\section triaxial-survey Survey of triaxial geodesics

Before delving into the nature of geodesics on a triaxial geodesic, it
is worth reviewing geodesics on an ellipsoid of revolution.  There are
two classes of simple closed geodesics (i.e., geodesics which close on
themselves without intersection): the equator and all the meridians.
All other geodesics oscillate between two equal and opposite circles of
latitude; but after completing a full oscillation in latitude these fall
slightly short (for an oblate ellipsoid) of completing a full circuit in
longitude.

Turning to the triaxial case, we find that there are only 3 simple
closed geodesics, the three principal sections of the ellipsoid given by
\f$x = 0\f$, \f$y = 0\f$, and \f$z = 0\f$.  To survey the other
geodesics, it is convenient to consider geodesics which intersect the
middle principal section, \f$y = 0\f$, at right angles.  Such geodesics
are shown in Figs. 2--6, where I use the same ellipsoid parameters as in
Fig. 1 and the same viewing direction.  In addition, The three principal
ellipses are shown in red in each of these figures.

If the starting point is \f$\beta_1 \in (-90^\circ, 90^\circ)\f$,
\f$\omega_1 = 0\f$, and \f$\alpha_1 = 90^\circ\f$, then the geodesic
encircles the ellipsoid in a "circumpolar" sense.  The geodesic
oscillates north and south of the equator; on each oscillation it
completes slightly less that a full circuit around the ellipsoid
resulting in the geodesic filling the area bounded by the two latitude
lines \f$\beta = \pm \beta_1\f$.  Two examples are given in
Figs. 2 and 3.  Figure 2 shows practically the same behavior as for an
oblate ellipsoid of revolution (because \f$a \approx b\f$).  However, if
the starting point is at a higher latitude (Fig. 3) the distortions
resulting from \f$a \ne b\f$ are evident.

<center>
<img src="triaxial-equatorial-a.png"
width=419 height=397 alt="Example of a circumpolar geodesic on a
triaxial ellipsoid">
Fig. 2
</center>\n
Fig. 2:
Example of a circumpolar geodesic on a triaxial ellipsoid.  The starting
point of this geodesic is \f$\beta_1 = 45.1^\circ\f$, \f$\omega_1 =
0^\circ\f$, and \f$\alpha_1 = 90^\circ\f$.

<center>
<img src="triaxial-equatorial-b.png"
width=419 height=397 alt="Another example of a circumpolar geodesic on a
triaxial ellipsoid">
Fig. 3
</center>\n
Fig. 3:
Another example of a circumpolar geodesic on a triaxial ellipsoid.  The
starting point of this geodesic is \f$\beta_1 = 87.48^\circ\f$, \f$\omega_1 =
0^\circ\f$, and \f$\alpha_1 = 90^\circ\f$.

If the starting point is \f$\beta_1 = 90^\circ\f$, \f$\omega_1 \in
(0^\circ, 180^\circ)\f$, and \f$\alpha_1 = 0^\circ\f$, then the geodesic
encircles the ellipsoid in a "transpolar" sense.  The geodesic
oscillates east and west of the ellipse \f$x = 0\f$; on each oscillation
it completes slightly more that a full circuit around the ellipsoid
resulting in the geodesic filling the area bounded by the two longitude
lines \f$\omega = \omega_1\f$ and \f$\omega = 180^\circ - \omega_1\f$.
If \f$a = b\f$, all meridians are geodesics; the effect of \f$a \ne b\f$
causes such geodesics to oscillate east and west.  Two examples are
given in Figs. 4 and 5.

<center>
<img src="triaxial-polar-a.png"
width=419 height=397 alt="Example of a transpolar geodesic on a
triaxial ellipsoid">
Fig. 4
</center>\n
Fig. 4:
Example of a transpolar geodesic on a triaxial ellipsoid.  The
starting point of this geodesic is \f$\beta_1 = 90^\circ\f$, \f$\omega_1 =
39.9^\circ\f$, and \f$\alpha_1 = 0^\circ\f$.

<center>
<img src="triaxial-polar-b.png"
width=419 height=397 alt="Another example of a transpolar geodesic on a
triaxial ellipsoid">
Fig. 5
</center>\n
Fig. 5:
Another example of a transpolar geodesic on a triaxial ellipsoid.  The
starting point of this geodesic is \f$\beta_1 = 90^\circ\f$, \f$\omega_1 =
9.966^\circ\f$, and \f$\alpha_1 = 0^\circ\f$.

If the starting point is \f$\beta_1 = 90^\circ\f$, \f$\omega_1 =
0^\circ\f$ (an umbilical point), and \f$\alpha_1 = 45^\circ\f$ (the
geodesic leaves the ellipse \f$y = 0\f$ at right angles), then the
geodesic repeatedly intersects the opposite umbilical point and returns to
its starting point.  However on each circuit the angle at which
intersects \f$y = 0\f$ becomes closer to \f$0^\circ\f$ or
\f$180^\circ\f$ so that asymptotically the geodesic lies on the ellipse
\f$y = 0\f$.  This is shown in Fig. 6.  Note that a single geodesic does
not fill an area on the ellipsoid.

<center>
<img src="triaxial-umbilical.png"
width=419 height=397 alt="Example of an umbilical geodesic on a
triaxial ellipsoid">
Fig. 6
</center>\n
Fig. 6:
Example of an umbilical geodesic on a triaxial ellipsoid.  The
starting point of this geodesic is \f$\beta_1 = 90^\circ\f$, \f$\omega_1 =
0^\circ\f$, and \f$\alpha_1 = 45^\circ\f$ and the geodesics is followed
forwards and backwards until it lies close to the plane \f$y = 0\f$ in
both directions.

Umbilical geodesic enjoy several interesting properties.
 - Through any point on the ellipsoid, there are two umbilical geodesics.
 - The geodesic distance between opposite umbilical points is the same
   regardless of the initial direction of the geodesic.
 - Whereas the closed geodesics on the ellipses \f$x = 0\f$ and \f$z =
   0\f$ are stable (an geodesic initially close to and nearly parallel to
   the ellipse remains close to the ellipse), the closed geodesic on the
   ellipse \f$y = 0\f$, which goes through all 4 umbilical points, is \e
   unstable.  If it is perturbed, it will swing out of the plane \f$y =
   0\f$ and flip around before returning to close to the plane.  (This
   behavior may repeat depending on the nature of the initial
   perturbation.).

\section triaxial-stab The stability of closed geodesics

The stability of the three simple closed geodesics can be determined by
examining the properties of Jacobi's solution.  In particular the
unstable behavior of umbilical geodesics was shown by Hart (1849).
However an alternative approach is to use the equations that Gauss
(1828) gives for a perturbed geodesic
\f[
\frac {d^2m}{ds^2} + Km = 0
\f]
where \f$m\f$ is the distance of perturbed geodesic from a reference
geodesic and \f$K\f$ is the Gaussian curvature of the surface.  If the
reference geodesic is closed, then this is a linear homogeneous
differential equation with periodic coefficients.  In fact it's a
special case of Hill's equation which can be treated using Floquet
theory, see <a href="http://dlmf.nist.gov/28.29">DLMF, &sect;28.29</a>.
Using the notation of &sect;3 of
<a href="http://dx.doi.org/10.1007/s00190-012-0578-z"> Algorithms for
geodesics</a>, the stability is determined by computing the reduced
length \f$m_{12}\f$ and the geodesic scales \f$M_{12}, M_{21}\f$ over
half the perimeter of the ellipse and determining the eigenvalues
\f$\lambda_{1,2}\f$ of
\f[
{\cal M} = \left(\begin{array}{cc}
M_{12} & m_{12}\\
-\frac{1 - M_{12}M_{21}}{m_{12}} & M_{21}
\end{array}\right).
\f]
Because \f$\mathrm{det}\,{\cal M} = 1\f$, the eigenvalues are determined
by \f$\mathrm{tr}\,{\cal M}\f$.  In particular if
\f$\left|\mathrm{tr}\,{\cal M}\right| < 2\f$, we have
\f$\left|\lambda_{1,2}\right| = 1\f$ and the solution is stable; if
\f$\left|\mathrm{tr}\,{\cal M}\right| > 2\f$, one of
\f$\left|\lambda_{1,2}\right|\f$ is larger than unity and the solution
is (exponentially) unstable.  In the transition case,
\f$\left|\mathrm{tr}\,{\cal M}\right| = 2\f$, the solution is stable
provided that the off-diagonal elements of \f${\cal M}\f$ are zero;
otherwise the solution is linearly unstable.

The exponential instability of the geodesic on the ellipse \f$y = 0\f$
is confirmed by this analysis and results from the resonance between the
natural frequency of the equation for \f$m\f$ and the driving frequency
when \f$b\f$ lies in \f$(c, a)\f$.  If \f$b\f$ is equal to either of the
other axes (and the triaxial ellipsoid degenerates to an ellipsoid of
revolution), then the solution is linearly unstable.  (For example, a
geodesic is which is close to a meridian on an oblate ellipsoid, slowly
moves away from that meridian.)

\section triaxial-inverse The inverse problem

In order to solve the inverse geodesic problem, it helps to have an
understanding of the properties of all the geodesics emanating from a
single point \f$(\beta_1, \omega_1)\f$.
 - If the point is an umbilical point, all the lines meet at the
   opposite umbilical point.
 - Otherwise, the first envelope of the geodesics is a 4-pointed
   astroid.  The cusps of the astroid lie on either \f$\beta = -
   \beta_1\f$ or \f$\omega = \omega_1 + \pi\f$; see
   <a href="http://dx.doi.org/10.1080/10586458.2003.10504515"> Sinclair
   (2003)</a>.
 - All geodesics intersect (or, in the case of \f$\alpha_1 = 0\f$ or
   \f$\pi\f$, touch) the line \f$\omega = \omega_1 + \pi\f$.
 - All geodesics intersect (or, in the case of \f$\alpha_1 =
   \pm\pi/2\f$, touch) the line \f$\beta = -\beta_1\f$.
 - Two geodesics with azimuths \f$\pm\alpha_1\f$ first intersect on
   \f$\omega = \omega_1 + \pi\f$ and their lengths to the point of
   intersection are equal.
 - Two geodesics with azimuths \f$\alpha_1\f$ and \f$\pi-\alpha_1\f$
   first intersect on \f$\beta = -\beta_1\f$ and their lengths to the
   point of intersection are equal.
 .
(These assertions follow directly from the equations for the geodesics;
some of them are somewhat surprising given the asymmetries of the
ellipsoid.)  Consider now terminating the geodesics from \f$(\beta_1,
\omega_1)\f$ at the point where they first intersect (or touch) the line
\f$\beta = -\beta_1\f$.  To focus the discussion, take \f$\beta_1 \le
0\f$.
 - The geodesics completely fill the portion of the ellipsoid satisfying
   \f$\beta \le -\beta_1\f$.
 - None of geodesics intersect any other geodesics.
 - Any initial portion of these geodesics is a shortest path.
 - Each geodesic intersects the line \f$\beta = \beta_2\f$, where
   \f$\beta_1 < \beta_2 < -\beta_1\f$, exactly once.
 - For a given \f$\beta_2\f$, this defines a continuous monotonic
   mapping of the circle of azimuths \f$\alpha_1\f$ to the circle of
   longitudes \f$\omega_2\f$.
 - If \f$\beta_2 = \pm \beta_1\f$, then the previous two assertions need
   to be modified similarly to the case for an ellipsoid of revolution.

These properties show that the inverse problem can be solved using
techniques similar to those employed for an ellipsoid of revolution (see
&sect;4 of
<a href="http://dx.doi.org/10.1007/s00190-012-0578-z"> Algorithms for
geodesics</a>).
 - If the points are opposite umbilical points, an arbitrary
   \f$\alpha_1\f$ may be chosen.
 - If the points are neighboring umbilical points, the shortest path
   lies on the ellipse \f$y = 0\f$.
 - If only one point is an umbilicial point, the azimuth at the
   non-umbilical point is found using the generalization of Clairaut's
   equation (given above) with \f$\gamma = 0\f$.
 - If both points lie on the equator \f$\beta = 0\f$, then determine the
   reduced length \f$m_{12}\f$ for the geodesic which is the shorter
   path along the ellipse \f$z = 0\f$.  If \f$m_{12} \ge 0\f$, then this
   is the shortest path on the ellipsoid; otherwise proceed to the
   general case (next).
 - Swap the points, if necessary, so that the first point is the one
   closest to a pole.  Estimate \f$\alpha_1\f$ (by some means) and solve
   the \e hybrid problem, i.e., determine the longitude \f$\omega_2\f$
   corresponding to the first intersection of the geodesic with \f$\beta
   = \beta_2\f$.  Adjust \f$\alpha_1\f$ so that the value of
   \f$\omega_2\f$ matches the given \f$\omega_2\f$ (there is a single
   root).  If a sufficiently close solution can be found, Newton's
   method can be employed since the necessary derivative can be
   expressed in terms of the reduced length \f$m_{12}\f$.

The shortest path found by this method is unique unless:
 - The length of the geodesic vanishes \f$s_{12}=0\f$, in which case any
   constant can be added to the azimuths.
 - The points are opposite umbilical points.  In this case,
   \f$\alpha_1\f$ can take on any value and \f$\alpha_2\f$ needs to be
   adjusted to maintain the value of \f$\tan\alpha_1 / \tan\alpha_2\f$.
 - \f$\beta_1 + \beta_2 = 0\f$ and \f$\cos\alpha_1\f$ and
   \f$\cos\alpha_2\f$ have opposite signs.  In this case, there another
   shortest geodesic with azimuths \f$\pi - \alpha_1\f$ and
   \f$\pi - \alpha_2\f$.

<center>
Back to \ref geodesic.  Forward to \ref transversemercator.  Up to
\ref contents.
</center>
**********************************************************************/
/**
\page transversemercator Transverse Mercator projection

<center>
Back to \ref triaxial.  Forward to \ref geocentric.  Up to \ref contents.
</center>

GeographicLib::TransverseMercator and
GeographicLib::TransverseMercatorExact provide accurate implementations
of the transverse Mercator projection.  The
 <a href="TransverseMercatorProj.1.html">TransverseMercatorProj</a>
utility provides an interface to these classes.

Go to
 - \ref testmerc
 - \ref tmseries
 - \ref tmfigures

References
 - L. Kr&uuml;ger,
   <a href="http://dx.doi.org/10.2312/GFZ.b103-krueger28"> Konforme
   Abbildung des Erdellipsoids in der Ebene</a> (Conformal mapping of
   the ellipsoidal earth to the plane), Royal Prussian Geodetic Institute,
   New Series 52, 172 pp. (1912).
 - L. P. Lee,
   Conformal Projections Based on Elliptic Functions,
   (B. V. Gutsell, Toronto, 1976), 128pp.,
   ISBN: 0919870163
   (Also appeared as:
   Monograph 16, Suppl. No. 1 to Canadian Cartographer, Vol 13).
   Part V, pp. 67--101,
   <a href="http://dx.doi.org/10.3138/X687-1574-4325-WM62">Conformal
   Projections Based On Jacobian Elliptic Functions</a>.
 - C. F. F. Karney,
   <a href="http://dx.doi.org/10.1007/s00190-011-0445-3">
   Transverse Mercator with an accuracy of a few nanometers,</a>
   J. Geodesy 85(8), 475--485 (Aug. 2011);
   preprint:
   <a href="http://arxiv.org/abs/1002.1417"> arXiv:1002.1417</a>;
   resource page:
   <a href="http://geographiclib.sf.net/tm.html"> tm.html</a>.
 .
The algorithm for GeographicLib::TransverseMercator is based on
Kr&uuml;ger (1912); that for GeographicLib::TransverseMercatorExact is
based on Lee (1976).

See <a href="http://geographiclib.sourceforge.net/tm-grid.kmz"
type="application/vnd.google-earth.kmz"> tm-grid.kmz</a>, for an
illustration of the exact transverse Mercator grid in Google Earth.

\section testmerc Test data for the transverse Mercator projection

A test set for the transverse Mercator projection is available at
 - <a href="http://sf.net/projects/geographiclib/files/testdata/TMcoords.dat.gz">
   TMcoords.dat.gz</a>
 .
This is about 17 MB (compressed).  This test set consists of a set of
geographic coordinates together with the corresponding transverse
Mercator coordinates.  The WGS84 ellipsoid is used, with central
meridian 0&deg;, central scale factor 0.9996 (the UTM value),
false easting = false northing = 0 m.

Each line of the test set gives 6 space delimited numbers
 - latitude (degrees, exact)
 - longitude (degrees, exact &mdash; see below)
 - easting (meters, accurate to 0.1 pm)
 - northing (meters, accurate to 0.1 pm)
 - meridian convergence (degrees, accurate to 10<sup>&minus;18</sup> deg)
 - scale (accurate to 10<sup>&minus;20</sup>)
 .
The latitude and longitude are all multiples of 10<sup>&minus;12</sup>
deg and should be regarded as exact, except that longitude =
82.63627282416406551 should be interpreted as exactly 90 (1 - \e e)
degrees.  These results are computed using Lee's formulas with
<a href="http://en.wikipedia.org/wiki/Maxima_(software)">Maxima</a>'s
bfloats and fpprec set to 80 (so the errors in the data are probably 1/2
of the values quoted above).  The Maxima code,
<a href="tm.mac">tm.mac</a> and <a href="ellint.mac">ellint.mac</a>,
used to prepare this data set is included in the distribution.  You will
need to have Maxima installed to use this code.  The comments at the top
of <a href="tm.mac">tm.mac</a> illustrate how to run it.

The contents of the file are as follows:
 - 250000 entries randomly distributed in lat in [0, 90], lon in [0, 90];
 - 1000 entries randomly distributed on lat in [0, 90], lon = 0;
 - 1000 entries randomly distributed on lat = 0, lon in [0, 90];
 - 1000 entries randomly distributed on lat in [0, 90], lon = 90;
 - 1000 entries close to lat = 90 with lon in [0, 90];
 - 1000 entries close to lat = 0, lon = 0 with lat &ge; 0, lon &ge; 0;
 - 1000 entries close to lat = 0, lon = 90 with lat &ge; 0, lon &le; 90;
 - 2000 entries close to lat = 0, lon = 90 (1 &minus; \e e) with lat &ge; 0;
 - 25000 entries randomly distributed in lat
   in [&minus;89, 0], lon in [90 (1 &minus; \e e), 90];
 - 1000 entries randomly distributed on lat in [&minus;89, 0], lon = 90;
 - 1000 entries randomly distributed on lat in [&minus;89, 0], lon = 90
   (1 &minus; \e e);
 - 1000 entries close to lat = 0, lon = 90 (lat < 0, lon &le; 90);
 - 1000 entries close to lat = 0, lon = 90 (1 &minus; \e e)
   (lat < 0, lon &le; 90 (1 &minus; \e e));
 .
(a total of 287000 entries).  The entries for lat < 0&deg; and
lon in [90&deg; (1 &minus; \e e), 90&deg;] use the "extended"
domain for the transverse Mercator projection explained in Sec. 5 of
<a href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>.  The first
258000 entries have lat &ge; 0&deg; and are suitable for testing
implementations following the standard convention.

\section tmseries Series approximation for transverse Mercator

Kr&uuml;ger (1912) gives a 4th-order approximation to the transverse
Mercator projection.  This is accurate to about 200 nm (200
nanometers) within the UTM domain.  Here we present the series
extended to 10th order.  By default, GeographicLib::TransverseMercator
uses the 6th-order approximation. The preprocessor macro
GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER can be used to select an order
from 4 thru 8.  The series expanded to order <i>n</i><sup>30</sup> are
given in <a href="tmseries30.html"> tmseries30.html</a>.

In the formulas below ^ indicates exponentiation (\e n^3 = \e n*\e n*\e
n) and / indicates real division (3/5 = 0.6).  The equations need to be
converted to Horner form, but are here left in expanded form so that
they can be easily truncated to lower order in \e n.  Some of the
integers here are not representable as 32-bit integers and will need to
be included as double literals.

\e A in Kr&uuml;ger, p. 12, eq. (5)
\verbatim
  A = a/(n + 1) * (1 + 1/4 * n^2
                     + 1/64 * n^4
                     + 1/256 * n^6
                     + 25/16384 * n^8
                     + 49/65536 * n^10);
\endverbatim

&gamma; in Kr&uuml;ger, p. 21, eq. (41)
\verbatim
alpha[1] =   1/2 * n
           - 2/3 * n^2
           + 5/16 * n^3
           + 41/180 * n^4
           - 127/288 * n^5
           + 7891/37800 * n^6
           + 72161/387072 * n^7
           - 18975107/50803200 * n^8
           + 60193001/290304000 * n^9
           + 134592031/1026432000 * n^10;
alpha[2] =   13/48 * n^2
           - 3/5 * n^3
           + 557/1440 * n^4
           + 281/630 * n^5
           - 1983433/1935360 * n^6
           + 13769/28800 * n^7
           + 148003883/174182400 * n^8
           - 705286231/465696000 * n^9
           + 1703267974087/3218890752000 * n^10;
alpha[3] =   61/240 * n^3
           - 103/140 * n^4
           + 15061/26880 * n^5
           + 167603/181440 * n^6
           - 67102379/29030400 * n^7
           + 79682431/79833600 * n^8
           + 6304945039/2128896000 * n^9
           - 6601904925257/1307674368000 * n^10;
alpha[4] =   49561/161280 * n^4
           - 179/168 * n^5
           + 6601661/7257600 * n^6
           + 97445/49896 * n^7
           - 40176129013/7664025600 * n^8
           + 138471097/66528000 * n^9
           + 48087451385201/5230697472000 * n^10;
alpha[5] =   34729/80640 * n^5
           - 3418889/1995840 * n^6
           + 14644087/9123840 * n^7
           + 2605413599/622702080 * n^8
           - 31015475399/2583060480 * n^9
           + 5820486440369/1307674368000 * n^10;
alpha[6] =   212378941/319334400 * n^6
           - 30705481/10378368 * n^7
           + 175214326799/58118860800 * n^8
           + 870492877/96096000 * n^9
           - 1328004581729009/47823519744000 * n^10;
alpha[7] =   1522256789/1383782400 * n^7
           - 16759934899/3113510400 * n^8
           + 1315149374443/221405184000 * n^9
           + 71809987837451/3629463552000 * n^10;
alpha[8] =   1424729850961/743921418240 * n^8
           - 256783708069/25204608000 * n^9
           + 2468749292989891/203249958912000 * n^10;
alpha[9] =   21091646195357/6080126976000 * n^9
           - 67196182138355857/3379030566912000 * n^10;
alpha[10]=   77911515623232821/12014330904576000 * n^10;
\endverbatim

&beta; in Kr&uuml;ger, p. 18, eq. (26*)
\verbatim
 beta[1] =   1/2 * n
           - 2/3 * n^2
           + 37/96 * n^3
           - 1/360 * n^4
           - 81/512 * n^5
           + 96199/604800 * n^6
           - 5406467/38707200 * n^7
           + 7944359/67737600 * n^8
           - 7378753979/97542144000 * n^9
           + 25123531261/804722688000 * n^10;
 beta[2] =   1/48 * n^2
           + 1/15 * n^3
           - 437/1440 * n^4
           + 46/105 * n^5
           - 1118711/3870720 * n^6
           + 51841/1209600 * n^7
           + 24749483/348364800 * n^8
           - 115295683/1397088000 * n^9
           + 5487737251099/51502252032000 * n^10;
 beta[3] =   17/480 * n^3
           - 37/840 * n^4
           - 209/4480 * n^5
           + 5569/90720 * n^6
           + 9261899/58060800 * n^7
           - 6457463/17740800 * n^8
           + 2473691167/9289728000 * n^9
           - 852549456029/20922789888000 * n^10;
 beta[4] =   4397/161280 * n^4
           - 11/504 * n^5
           - 830251/7257600 * n^6
           + 466511/2494800 * n^7
           + 324154477/7664025600 * n^8
           - 937932223/3891888000 * n^9
           - 89112264211/5230697472000 * n^10;
 beta[5] =   4583/161280 * n^5
           - 108847/3991680 * n^6
           - 8005831/63866880 * n^7
           + 22894433/124540416 * n^8
           + 112731569449/557941063680 * n^9
           - 5391039814733/10461394944000 * n^10;
 beta[6] =   20648693/638668800 * n^6
           - 16363163/518918400 * n^7
           - 2204645983/12915302400 * n^8
           + 4543317553/18162144000 * n^9
           + 54894890298749/167382319104000 * n^10;
 beta[7] =   219941297/5535129600 * n^7
           - 497323811/12454041600 * n^8
           - 79431132943/332107776000 * n^9
           + 4346429528407/12703122432000 * n^10;
 beta[8] =   191773887257/3719607091200 * n^8
           - 17822319343/336825216000 * n^9
           - 497155444501631/1422749712384000 * n^10;
 beta[9] =   11025641854267/158083301376000 * n^9
           - 492293158444691/6758061133824000 * n^10;
 beta[10]=   7028504530429621/72085985427456000 * n^10;
\endverbatim

The high-order expansions for &alpha; and &beta; were produced by the
Maxima program <a href="tmseries.mac"> tmseries.mac</a> (included in the
distribution).  To run, start Maxima and enter
\verbatim
  load("tmseries.mac")$
\endverbatim
Further instructions are included at the top of the file.

\section tmfigures Figures from paper on transverse Mercator projection

This section gives color versions of the figures in
 - C. F. F. Karney,
   <a href="http://dx.doi.org/10.1007/s00190-011-0445-3">
   Transverse Mercator with an accuracy of a few nanometers,</a>
   J. Geodesy 85(8), 475--485 (Aug. 2011);
   preprint:
   <a href="http://arxiv.org/abs/1002.1417"> arXiv:1002.1417</a>;
   resource page:
   <a href="http://geographiclib.sf.net/tm.html"> tm.html</a>.

\b NOTE:
The numbering of the figures matches that in the paper cited above.  The
figures in this section are relatively small in order to allow them to
be displayed quickly.  Vector versions of these figures are available in
<a href="http://geographiclib.sourceforge.net/tm-figs.pdf"
type="application/pdf"> tm-figs.pdf</a>.  This allows you to magnify the
figures to show the details more clearly.

<center>
<img src="gauss-schreiber-graticule-a.png"
width=667 height=305 alt="Gauss-Schreiber transverse Mercator graticule
(scaled)">
Fig. 1(a)
</center>\n
Fig. 1(a):
The graticule for the spherical transverse Mercator projection.
The equator lies on \e y = 0.
Compare this with Lee, Fig. 1 (right), which shows the graticule for
half a sphere, but note that in his notation \e x and \e y have switched
meanings.  The graticule for the ellipsoid differs from that for a
sphere only in that the latitude lines have shifted slightly.  (The
conformal transformation from an ellipsoid to a sphere merely relabels
the lines of latitude.)  This projection places the point latitude =
0&deg;, longitude = 90&deg; at infinity.

<center>
<img src="gauss-krueger-graticule-a.png"
width=667 height=305 alt="Gauss-Krueger transverse Mercator graticule
(scaled)">
Fig. 1(b)
</center>\n
Fig. 1(b):
The graticule for the Gauss-Kr&uuml;ger transverse Mercator projection.
The equator lies on \e y = 0 for longitude < 81&deg;; beyond
this, it arcs up to meet \e y = 1.  Compare this with Lee, Fig. 45
(upper), which shows the graticule for the International Ellipsoid.  Lee,
Fig. 46, shows the graticule for the entire ellipsoid.  This projection
(like the Thompson projection) projects the ellipsoid to a finite area.

<center>
<img src="thompson-tm-graticule-a.png"
width=667 height=305 alt="Thompson transverse Mercator graticule
(scaled)">
Fig. 1(c)
</center>\n
Fig. 1(c):
The graticule for the Thompson transverse Mercator projection.  The
equator lies on \e y = 0 for longitude < 81&deg;; at longitude =
81&deg;, it turns by 120&deg; and heads for \e y = 1.
Compare this with Lee, Fig. 43, which shows the graticule for the
International Ellipsoid.  Lee, Fig. 44, shows the graticule for the
entire ellipsoid.  This projection (like the Gauss-Kr&uuml;ger
projection) projects the ellipsoid to a finite area.

<center>
<img src="gauss-krueger-error.png"
width=672 height=552 alt="Truncation error for Gauss-Krueger series">
Fig. 2
</center>\n
Fig. 2:
The truncation error for the series for the Gauss-Kr&uuml;ger transverse
Mercator projection.  The blue curves show the truncation error for the
order of the series \e J = 2 (top) thru \e J = 12 (bottom).  The red
curves show the combined truncation and round-off errors for
 - float and \e J = 4 (top)
 - double and \e J = 6 (middle)
 - long double and \e J = 8 (bottom)

<center>
<img src="gauss-krueger-graticule.png"
width=673 height=323 alt="Gauss-Krueger transverse Mercator extended
graticule">\n
Fig. 3(a)
</center>
Fig. 3(a):
The graticule for the extended domain.  The blue lines show
latitude and longitude at multiples of 10&deg;.  The green lines
show 1&deg; intervals for longitude in [80&deg;, 90&deg;] and latitude in
[&minus;5&deg;, 10&deg;].

<center>
<img src="gauss-krueger-convergence-scale.png"
width=673 height=323 alt="Gauss-Krueger transverse Mercator
convergence and scale">
Fig. 3(b)
</center>\n
Fig. 3(b):
The convergence and scale for the Gauss-Kr&uuml;ger
transverse Mercator projection in the extended domain.  The blue lines
emanating from the top left corner (the north pole) are lines of
constant convergence.  Convergence = 0&deg; is given by the
dog-legged line joining the points (0,1), (0,0), (1.71,0),
(1.71,&minus;&infin;).
Convergence = 90&deg; is given by the line y = 1.  The other
lines show multiples of 10&deg; between 0&deg; and
90&deg;.  The other blue, the green and the black lines show
scale = 1 thru 2 at intervals of 0.1, 2 thru 15 at intervals of 1, and
15 thru 35 at intervals of 5.  Multiples of 5 are shown in black,
multiples of 1 are shown in blue, and the rest are shown in green.
Scale = 1 is given by the line segment (0,0) to (0,1).  The red line
shows the equator between lon = 81&deg; and 90&deg;.  The
scale and convergence at the branch point are 1/\e e = 10 and
0&deg;, respectively.

<center>
<img src="thompson-tm-graticule.png"
width=661 height=321 alt="Thompson transverse Mercator extended graticule">
Fig. 3(c)
</center>\n
Fig. 3(c):
The graticule for the Thompson transverse Mercator
projection for the extended domain.  The range of the projection is the
rectangular region shown
 - 0 &le; \e u &le; K(<i>e</i><sup>2</sup>),
   0 &le; \e v &le; K(1 &minus; <i>e</i><sup>2</sup>)
 .
The coloring of the lines is the same as Fig. 3(a), except that latitude
lines extended down to &minus;10&deg; and a red line has been added
showing the line \e y = 0 for \e x > 1.71 in the Gauss-Kr&uuml;ger
projection (Fig. 3(a)).  The extended Thompson projection figure has
reflection symmetry on all the four sides of Fig. 3(c).

<center>
Back to \ref triaxial.  Forward to \ref geocentric.  Up to \ref contents.
</center>
**********************************************************************/
/**
\page geocentric Geocentric coordinates

<center>
Back to \ref transversemercator.  Forward to \ref old.  Up to \ref contents.
</center>

The implementation of GeographicLib::Geocentric::Reverse is adapted from
 - H. Vermeille,
   <a href="http://dx.doi.org/10.1007/s00190-002-0273-6">
   Direct transformation from geocentric coordinates to geodetic
   coordinates</a>, J. Geodesy 76, 451--454 (2002).

This provides a closed-form solution but can't directly be applied close to
the center of the earth.  Several changes have been made to remove this
restriction and to improve the numerical accuracy.  Now the method is
accurate for all inputs (even if \e h is infinite).  The changes are
described in Appendix B of
 - C. F. F. Karney,
   <a href="http://arxiv.org/abs/1102.1215v1">Geodesics
   on an ellipsoid of revolution</a>,
   Feb. 2011; preprint
   <a href="http://arxiv.org/abs/1102.1215v1">arxiv:1102.1215v1</a>.

The problems encountered near the center of the ellipsoid are:
 - There's a potential division by zero in the definition of \e s. The
   equations are easily reformulated to avoid this problem.
 - <i>t</i><sup>3</sup> may be negative.  This is OK; we just take the
   real root.
 - The solution for \e t may be complex.  However this leads to 3 real roots
   for <i>u</i>/\e r.  It's then just a matter of picking the one that computes
   the geodetic result which minimizes |\e h| and which avoids large
   round-off errors.
 - Some of the equations result in a large loss of accuracy due to
   subtracting nearly equal quantities.  E.g., \e k= sqrt(\e u + \e v +
   <i>w</i><sup>2</sup>) &minus; \e w is inaccurate if \e u + \e v is small;
   we can fix this by writing \e k = (\e u + \e v)/(sqrt(\e u + \e v +
   <i>w</i><sup>2</sup>) + \e w).

The error is computed as follows.  Write a version of
Geocentric::WGS84.Forward which uses long doubles (including using long
doubles for the WGS84 parameters).  Generate random (long double)
geodetic coordinates (\e lat0, \e lon0, \e h0) and use the "long double"
WGS84.Forward to obtain the corresponding (long double) geocentric
coordinates (\e x0, \e y0, \e z0).  [We restrict \e h0 so that
\e h0 &ge; &minus; \e a (1 &minus; <i>e</i><sup>2</sup>) / sqrt(1 &minus;
<i>e</i><sup>2</sup> sin<sup>2</sup>\e lat0), which ensures that (\e
lat0, \e lon0, \e h0) is the principal geodetic inverse of (\e x0, \e
y0, \e z0).]  Because the forward calculation is numerically stable and
because long doubles (on Linux systems using g++) provide 11 bits
additional accuracy (about 3.3 decimal digits), we regard this set of
test data as exact.

Apply the double version of WGS84.Reverse to (\e x0, \e y0, \e z0) to
compute the approximate geodetic coordinates (\e lat1, \e lon1, \e h1).
Convert (\e lat1 &minus; \e lat0, \e lon1 &minus; \e lon0) to a
distance, \e ds, on the surface of the ellipsoid and define \e err =
hypot(\e ds, \e h1 &minus; \e h0).  For |\e h0| < 5000 km, we have \e
err < 7 nm (7 nanometers).

This methodology is not very useful very far from the globe, because the
absolute errors in the approximate geodetic height become large, or
within 50 km of the center of the earth, because of errors in computing
the approximate geodetic latitude.  To illustrate the second issue, the
maximum value of \e err for \e h0 < 0 is about 80 mm.  The error is
maximum close to the circle given by geocentric coordinates satisfying
hypot(\e x, \e y) = \e a <i>e</i><sup>2</sup> (= 42.7 km), \e z = 0.
(This is the center of meridional curvature for \e lat = 0.)  The
geodetic latitude for these points is \e lat = 0.  However, if we move 1
nm towards the center of the earth, the geodetic latitude becomes 0.04",
a distance of 1.4 m from the equator.  If, instead, we move 1 nm up, the
geodetic latitude becomes 7.45", a distance of 229 m from the equator.
In light of this, Reverse does quite well in this vicinity.

To obtain a practical measure of the error for the general case we define
- <i>err</i><sub>h</sub> = |\e h1 &minus; \e h0| / max(1, \e h0 / \e a)
- for \e h0 > 0, <i>err</i><sub>out</sub> = \e ds
- for \e h0 < 0, apply the long double version of WGS84.Forward to (\e
  lat1, \e lon1, \e h1) to give (\e x1, \e y1, \e z1) and compute
  <i>err</i><sub>in</sub> = hypot(\e x1 &minus; \e x0, \e y1 &minus; \e
  y0, \e z1 &minus; \e z0).
.
We then find <i>err</i><sub>h</sub> < 8 nm, <i>err</i><sub>out</sub> < 4 nm,
and <i>err</i><sub>in</sub> < 7 nm.  (1 nm = 1 nanometer.)

The testing has been confined to the WGS84 ellipsoid.  The method will work
for all ellipsoids used in terrestrial geodesy.  However, the central region,
which leads to multiple real roots for the cubic equation in Reverse, pokes
outside the ellipsoid (at the poles) for ellipsoids with \e e > 1/sqrt(2).
Reverse has not been analyzed for this case.  Similarly ellipsoids which are
very nearly spherical near yield inaccurate results due to underflow; in the
other hand, the case of the sphere, \e f = 0, is treated specially and gives
accurate results.

Other comparable methods are K. M. Borkowski,
<a href="http://dx.doi.org/10.1007/BF00643807"> Transformation
of geocentric to geodetic coordinates without approximations</a>,
Astrophys. Space Sci. 139, 1--4 (1987)
(<a href="http://dx.doi.org/10.1007/BF00656995"> erratum</a>)
and T. Fukushima,
<a href="http://dx.doi.org/10.1007/s001900050271"> Fast transform from
geocentric to geodetic coordinates</a>, J. Geodesy 73, 603--610 (1999).
However the choice of independent variables in these methods leads
to a loss of accuracy for points near the equatorial plane.

<center>
Back to \ref transversemercator.  Forward to \ref old.  Up to \ref contents.

</center>
**********************************************************************/
/**
\page old Old versions

<center>
Back to \ref geocentric.  Up to \ref contents.
</center>

List of versions in reverse chronological order together with a brief
list of changes.  (Note: Old versions of the library use a year-month
style of numbering.  Now, the library uses a major and minor version
number.)  Recent versions of %GeographicLib are available at
<a href="http://sf.net/projects/geographiclib/files/distrib/">
http://sourceforge.net/projects/geographiclib/files/distrib/</a>.
Older versions are in
<a href="http://sf.net/projects/geographiclib/files/distrib/archive/">
http://sourceforge.net/projects/geographiclib/files/distrib/archive/</a>.

The corresponding documentation for these versions is obtained by
clicking on the &ldquo;Version <i>m.nn</i>&rdquo; links below.  Some of
the links in the documentaion of older versions may be out of date (in
particular the links for the source code will not work if the code has
been migrated to the archive subdirectory).  All the releases are
available as tags &ldquo;r<i>m.nn</i>&rdquo; in the the "release" branch
of the git repository for %GeographicLib.

 - <a href="http://geographiclib.sf.net/1.30">Version 1.30</a>
   (released 2013-02-27)
   - Changes to geodesic routines:
     - fix BUG in fail-safe mechanisms in GeographicLib::Geodesic::Inverse;
     - the command line utility Geod is now called
       <a href="GeodSolve.1.html">GeodSolve</a>;
     - allow addition of polygon edges in GeographicLib::PolygonArea;
     - add full Maxima implementation of geodesic algorithms.

 - <a href="http://geographiclib.sf.net/1.29">Version 1.29</a>
   (released 2013-01-16)
   - Changes to allow compilation with libc++ (courtesy of Kal Conley).
   - Add description of geodesics on triaxial ellipsoid to
     documentation.
   - Update journal reference for "Algorithms for geodesics".

 - <a href="http://geographiclib.sf.net/1.28">Version 1.28</a>
   (released 2012-12-11)
   - Changes to geodesic routines:
     - compute longitude difference exactly;
     - hence fix BUG in area calculations for polygons with vertices very
       close to the prime meridian;
     - fix BUG is geoddistance.m where the value of m12 was wrong for
       meridional geodesics;
     - add Matlab implementations of the geodesic projections;
     - remove unneeded special code for geodesics which start at a pole;
     - include polygon area routine in C and Fortran implementations;
     - add doxygen documentation for C and Fortran libraries.

 - <a href="http://geographiclib.sf.net/1.27">Version 1.27</a>
   (released 2012-11-29)
   - Changes to geodesic routines:
     - add native Matlab implementations: geoddistance.m, geodreckon.m,
       geodarea.m;
     - add C and Fortran implementations;
     - improve the solution of the direct problem so that the series
       solution is accurate to round off for |<i>f</i>| &lt; 1/50;
     - tighten up the convergence criteria for solution of the inverse
       problem;
     - no longer signal failures of convergence with NaNs (a slightly
       less accurate answer is returned instead).
   - Fix GeographicLib::DMS::Decode double rounding BUG.
   - On MacOSX platforms with the cmake configuration, universal
     binaries are built.

 - <a href="http://geographiclib.sf.net/1.26">Version 1.26</a>
   (released 2012-10-22)
   - Replace the series used for geodesic areas by one with better
     convergence (this only makes an appreciable difference if
     |<i>f</i>| &gt; 1/150).

 - <a href="http://geographiclib.sf.net/1.25">Version 1.25</a>
   (released 2012-10-16)
   - Changes to geodesic calculations:
     - restart Newton's method in Geodesic::Inverse when it goes awry;
     - back up Newton's method with the bisection method;
     - GeographicLib::Geodesic::Inverse now converges for any value of \e f;
     - add GeographicLib::GeodesicExact and
       GeographicLib::GeodesicLineExact which are formulated in terms
       of elliptic integrals and thus yield accurate results even for
       very eccentric ellipsoids.
     - the -E option to <a href="GeodSolve.1.html">Geod</a> invokes these
       exact classes.
   - Add functionality to GeographicLib::EllipticFunction:
     - add all the traditional elliptic integrals;
     - remove restrictions on argument range for incomplete elliptic integrals;
     - allow imaginary modulus for elliptic integrals and elliptic functions;
     - make interface to the symmetric elliptic integrals public.
   - Allow GeographicLib::Ellipsoid to be copied.
   - Changes to the build tools:
     - cmake uses folders in Visual Studio to reduce clutter;
     - allow precision of reals to be set in cmake;
     - fail gracefully in the absence of pod documentation tools;
     - remove support for maintainer tasks in Makefile.mk.

 - <a href="http://geographiclib.sf.net/1.24">Version 1.24</a>
   (released 2012-09-22)
   - Allow the specification of the hemisphere in UTM coordinates in
     order to provide continuity across the equator:
     - add GeographicLib::UTMUPS::Transfer;
     - add GeographicLib::GeoCoords::UTMUPSRepresentation(bool, int) and
       GeographicLib::GeoCoords::AltUTMUPSRepresentation(bool, int);
     - use the hemisphere letter in, e.g.,
       <a href="GeoConvert.1.html">GeoConvert</a> -u -z 31N.
   - Add GeographicLib::UTMUPS::DecodeEPSG and
     GeographicLib::UTMUPS::EncodeEPSG.
   - cmake changes:
     - restore support for cmake 2.4.x;
     - explicitly check version of doxygen.
   - Fix building under cygwin.
   - Document restrictions on \e f in \ref intro.
   - Fix python interface to work with version 2.6.x.

 - <a href="http://geographiclib.sf.net/1.23">Version 1.23</a>
   (released 2012-07-17)
  - Documentation changes:
    - remove html documentation from distribution and use web links if
      doxygen is not available;
    - use doxygen tags to document exceptions;
    - begin migrating the documentation to using Greek letters where
      appropriate (requires doxygen 1.8.1.2 or later).
  - Add GeographicLib::Math::AngNormalize and
    GeographicLib::Math::AngNormalize2; the allowed range for longitudes
    and azimuths widened to [&minus;540&deg;, 540&deg;).
  - GeographicLib::DMS::Decode understands more unicode symbols.
  - GeographicLib::Geohash uses geohash code "nan" to stand for not a
    number.
  - Add GeographicLib::Ellipsoid::NormalCurvatureRadius.
  - Various fixes in GeographicLib::LambertConformalConic,
    GeographicLib::TransverseMercator,
    GeographicLib::PolarStereographic, and GeographicLib::Ellipsoid to
    handle reverse projections of points near infinity.
  - Fix programming blunder in GeographicLib::LambertConformalConic::Forward
    (incorrect results were returned if the tangent latitude was
    negative).

 - <a href="http://geographiclib.sf.net/1.22">Version 1.22</a>
   (released 2012-05-27)
  - Add GeographicLib::Geohash and GeographicLib::Ellipsoid classes.
  - Fix bug in GeographicLib::AlbersEqualArea for very prolate
    ellipsoids (<i>b</i><sup>2</sup> > 2 <i>a</i><sup>2</sup>).
  - cmake changes:
    - optionally use PACKAGE_PATH and INSTALL_PATH to determine
      CMAKE_INSTALL_PREFIX;
    - use COMMON_INSTALL_PATH to determine layout of installation
      directories;
    - as a consequence, the installation paths for the documentation,
      and python and matlab interfaces are shortened for Windows;
    - zip source distribution now uses DOS line endings;
    - the tests work in debug mode for Windows;
    - default setting of GEOGRAPHICLIB_DATA does not depend on
      CMAKE_INSTALL_PREFIX;
    - add a cmake configuration for build tree.

 - <a href="http://geographiclib.sf.net/1.21">Version 1.21</a>
   (released 2012-04-25)
  - Support colon-separated DMS output:
    - GeographicLib::DMS::Encode and
      GeographicLib::GeoCoords::DMSRepresentation generalized;
    - <a href="GeoConvert.1.html">GeoConvert</a> and
      <a href="GeodSolve.1.html">Geod</a> now accept a -: option.
  - <a href="GeoidEval.1.html">GeoidEval</a> does not print the gradient
    of the geoid height by default (because it's subject to large
    errors); give the -g option to get the gradient printed.
  - Work around optimization BUG in GeographicLib::Geodesic::Inverse
    with tdm mingw g++ version 4.6.1.
  - autoconf fixed to ensure that that out-of-sources builds work;
    document this as the preferred method of using autoconf.
  - cmake tweaks:
    - simplify the configuration of doxygen;
    - allow the Matlab compiler to be specified with the
      MATLAB_COMPILER option.

 - <a href="http://geographiclib.sf.net/1.20">Version 1.20</a>
   (released 2012-03-23)
  - cmake tweaks:
    - improve find_package's matching of compiler versions;
    - CMAKE_INSTALL_PREFIX set from CMAKE_PREFIX_PATH if available;
    - add "x64" to the package name for the 64-bit binary installer;
    - fix cmake warning with Visual Studio Express.
  - Fix GeographicLib::SphericalEngine to deal with aggessive iterator
    checking by Visual Studio.
  - Fix transcription BUG is Geodesic.js.

 - <a href="http://geographiclib.sf.net/1.19">Version 1.19</a>
   (released 2012-03-13)
  - Slight improvement in GeographicLib::Geodesic::Inverse for very
    short lines.
  - Fix argument checking tests in GeographicLib::MGRS::Forward.
  - Add --comment-delimiter and --line-separator options to the \ref
    utilities.
  - Add installer for 64-bit Windows; the compiled Matlab interface is
    supplied with the Windows 64-bit installer only.

 - <a href="http://geographiclib.sf.net/1.18">Version 1.18</a>
   (released 2012-02-18)
  - Improve documentation on configuration with cmake.
  - cmake's find_package ensures that the compiler versions match on Windows.
  - Improve documentation on compiling Matlab interface.
  - Binary installer for Windows installs under C:/pkg-vc10 by default.

 - <a href="http://geographiclib.sf.net/1.17">Version 1.17</a>
   (released 2012-01-21)
  - Work around optimization BUG in GeographicLib::Geodesic::Inverse
    with g++ version 4.4.0 (mingw).
  - Fix BUG in argument checking with GeographicLib::OSGB::GridReference.
  - Fix missing include file in GeographicLib::SphericalHarmonic2.
  - Add simple examples of usage for each class.
  - Add internal documentation to the cmake configuration files.

 - <a href="http://geographiclib.sf.net/1.16">Version 1.16</a>
   (released 2011-12-07)
  - Add calculation of the earth's gravitational field:
    - add GeographicLib::NormalGravity GeographicLib::GravityModel and
      GeographicLib::GravityCircle classes;
    - add command line utility
      <a href="Gravity.1.html">Gravity</a>;
    - add \ref gravity;
    - add GeographicLib::Constants::WGS84_GM(),
      GeographicLib::Constants::WGS84_omega(), and similarly for GRS80.
  - Build uses GEOGRAPHICLIB_DATA to specify a common parent directory
     for geoid, gravity, and magnetic data (instead of
     GEOGRAPHICLIB_GEOID_PATH, etc.); similarly,
     <a href="GeoidEval.1.html">GeoidEval</a>,
     <a href="Gravity.1.html">Gravity</a>, and
     <a href="MagneticField.1.html">MagneticField</a>, look at the
     environment variable GEOGRAPHICLIB_DATA to locate the data.
  - Spherical harmonic software changes:
    - capitalize enums GeographicLib::SphericalHarmonic::FULL and
      GeographicLib::SphericalHarmonic::SCHMIDT (the lower case names
      are retained but deprecated);
    - optimize the sum by using a static table of square roots which is
      updated by GeographicLib::SphericalEngine::RootTable;
    - avoid overflow for high degree models.
  - Magnetic software fixes:
    - fix documentation BUG in GeographicLib::MagneticModel::Circle;
    - make GeographicLib::MagneticModel constructor explicit;
    - provide default GeographicLib::MagneticCircle constructor;
    - add additional inspector functions to
      GeographicLib::MagneticCircle;
    - add -c option to <a href="MagneticField.1.html">MagneticField</a>;
    - default height to zero in
      <a href="MagneticField.1.html">MagneticField</a>.

 - <a href="http://geographiclib.sf.net/1.15">Version 1.15</a>
   (released 2011-11-08)
  - Add calculation of the earth's magnetic field:
    - add GeographicLib::MagneticModel and GeographicLib::MagneticCircle
      classes;
    - add command line utility
      <a href="MagneticField.1.html">MagneticField</a>;
    - add \ref magnetic;
    - add \ref magneticinst;
    - add \ref magneticformat;
    - add classes GeographicLib::SphericalEngine,
      GeographicLib::CircularEngine, GeographicLib::SphericalHarmonic,
      GeographicLib::SphericalHarmonic1, and
      GeographicLib::SphericalHarmonic2. which sum spherical harmonic
      series.
  - Add GeographicLib::Utility class to support I/O and date
    manipulation.
  - Cmake configuration includes a _d suffix on the library built in
    debug mode.
  - For the Python package, include manifest and readme files; don't
    install setup.py for non-Windows systems.
  - Include Doxygen tag file in distribution as doc/html/Geographic.tag.

 - <a href="http://geographiclib.sf.net/1.14">Version 1.14</a>
   (released 2011-09-30)
  - Ensure that geographiclib-config.cmake is relocatable.
  - Allow more unicode symbols to be used in GeographicLib::DMS::Decode.
  - Modify <a href="GeoidEval.1.html">GeoidEval</a> so that it can be
    used to convert the height datum for LIDAR data.
  - Modest speed-up of Geodesic::Inverse.
  - Changes in python interface:
    - FIX BUG in transcription of Geodesic::Inverse;
    - include setup.py for easy installation;
    - python only distribution is available at
      http://pypi.python.org/pypi/geographiclib
  - Supply a minimal Qt qmake project file for library
    src/Geographic.pro.

 - <a href="http://geographiclib.sf.net/1.13">Version 1.13</a>
   (released 2011-08-13)
  - Changes to I/O:
    - allow : (colon) to be used as a DMS separator in
      GeographicLib::DMS::Decode(const std::string&, flag&);
    - also accept Unicode symbols for degrees, minutes, and seconds
      (coded as UTF-8);
    - provide optional \e swaplatlong argument to various
      GeographicLib::DMS and GeographicLib::GeoCoords functions to make
      longitude precede latitude;
    - <a href="GeoConvert.1.html">GeoConvert</a> now has a -w option to
      make longitude precede latitude on input and output;
    - include a Javascript version of GeographicLib::DMS.
  - Slight improvement in starting guess for solution of geographic
    latitude in terms of conformal latitude in TransverseMercator,
    TransverseMercatorExact, and LambertConformalConic.
  - For most classes, get rid of const member variables so that the
    default copy assignment works.
  - Put GeographicLib::Math and GeographicLib::Accumulator in their own
    header files.
  - Remove unused "fast" GeographicLib::Accumulator method.
  - Reorganize the \ref python.
  - Withdraw some deprecated routines.
  - cmake changes:
    - include FindGeographic.cmake in distribution;
    - building with cmake creates and installs
      geographiclib-config.cmake;
    - better support for building a shared library under Windows.

 - <a href="http://geographiclib.sf.net/1.12">Version 1.12</a>
   (released 2011-07-21)
  - Change license to MIT/X11.
  - Add GeographicLib::PolygonArea class and equivalent Matlab function.
  - Provide Javascript and Python implementations of geodesic routines.
  - Fix Windows installer to include runtime dlls for Matlab.
  - Fix (innocuous) unassigned variable in Geodesic::GenInverse.
  - Geodesic routines in Matlab return a12 as first column of aux return
    value (incompatible change).
  - A couple of code changes to enable compilation with Visual
    Studio 2003.

 - <a href="http://geographiclib.sf.net/1.11">Version 1.11</a>
   (released 2011-06-27)
  - Changes to <a href="Planimeter.1.html">Planimeter</a>:
    - add -l flag to <a href="Planimeter.1.html">Planimeter</a> for polyline
      calculations;
    - trim precision of area to 3 decimal places;
    - FIX BUG with pole crossing edges (due to compiler optimization).
  - <a href="GeodSolve.1.html">Geod</a> no longer reports the reduced
    length by default; however the -f flag still reports this and in
    addition gives the geodesic scales and the geodesic area.
  - FIX BUGS (compiler-specific) in inverse geodesic calculations.
  - FIX BUG: accommodate tellg() returning &minus;1 at end of string.
  - Change way flattening of the ellipsoid is specified:
    - constructors take \e f argument which is taken to be the
      flattening if \e f &lt; 1 and the inverse flattening otherwise
      (this is a compatible change for spheres and oblate ellipsoids, but it
      is an INCOMPATIBLE change for prolate ellipsoids);
    - the -e arguments to the \ref utilities are handled similarly; in
      addition, simple fractions, e.g., 1/297, can be used for the
      flattening;
    - introduce GeographicLib::Constants::WGS84_f() for the WGS84
      flattening (and deprecate Constants::WGS84_r() for the inverse
      flattening);
    - most classes have a Flattening() member function;
    - InverseFlattening() has been deprecated (and now returns inf for a
      sphere, instead of 0).

 - <a href="http://geographiclib.sf.net/1.10">Version 1.10</a>
   (released 2011-06-11)
  - Improvements to Matlab/Octave interface:
    - add {geocentric,localcartesian}{forward,reverse};
    - make geographiclibinterface more general;
    - install the source for the interface;
    - cmake compiles the interface if ENABLE_MATLAB=ON;
    - include compiled interface with Windows binary installer.
  - Fix various configuration issues
    - autoconf did not install Config.h;
    - cmake installed in man/man1 instead of share/man/man1;
    - cmake did not set the rpath on the tools.

 - <a href="http://geographiclib.sf.net/1.9">Version 1.9</a>
   (released 2011-05-28)
  - FIX BUG in area returned by
    <a href="Planimeter.1.html">Planimeter</a> for pole encircling polygons.
  - FIX BUG in error message reported when DMS::Decode reads the string
    "5d.".
  - FIX BUG in AlbersEqualArea::Reverse (lon0 not being used).
  - Ensure that all exceptions thrown in the \ref utilities are caught.
  - Avoid using catch within GeographicLib::DMS.
  - Move GeographicLib::Accumulator class from Planimeter.cpp to
    Constants.hpp.
  - Add GeographicLib::Math::sq<T>.
  - Simplify \ref geoidinst
    - add geographiclib-get-geoids for Unix-like systems;
    - add installers for Windows.
  - Provide cmake support:
    - build binary installer for Windows;
    - include regression tests;
    - add --input-string, --input-file, --output-file options to the
      \ref utilities to support tests.
  - Rename utility EquidistantTest as <a href="GeodesicProj.1.html">
    GeodesicProj</a> and TransverseMercatorTest as
    <a href="TransverseMercatorProj.1.html"> TransverseMercatorProj</a>.
  - Add <a href="ConicProj.1.html"> ConicProj</a>.
  - Reverse the initial sense of the -s option for
    <a href="Planimeter.1.html"> Planimeter</a>.
  - Migrate source from subversion to
    <a href="http://geographiclib.git.sf.net/git/gitweb-index.cgi">
    git</a>.

 - <a href="http://geographiclib.sf.net/1.8">Version 1.8</a>
   (released 2011-02-22)
  - Optionally return rotation matrix from GeographicLib::Geocentric and
    GeographicLib::LocalCartesian.
  - For the \ref utilities, supply man pages, -h prints the synopsis,
    --help prints the man page, --version prints the version.
  - Use accurate summation in <a href="Planimeter.1.html">Planimeter</a>.
  - Add 64-bit targets for Visual Studio 2010.
  - Use templates for defining math functions and some constants.
  - GeographicLib::Geoid updates
    - Add GeographicLib::Geoid::DefaultGeoidPath and
      GeographicLib::Geoid::DefaultGeoidName;
    - <a href="GeoidEval.1.html">GeoidEval</a> uses environment variable
      GEOID_NAME as the default geoid;
    - Add --msltohae and --haetomsl as
      <a href="GeoidEval.1.html">GeoidEval</a> options (and don't
      document the single hyphen versions).
  - Remove documentation that duplicates papers on transverse Mercator
    and geodesics.

 - <a href="http://geographiclib.sf.net/1.7">Version 1.7</a>
   (released 2010-12-21)
  - FIX BUG in scale returned by GeographicLib::LambertConformalConic::Reverse.
  - Add GeographicLib::AlbersEqualArea projection.
  - Library created by Visual Studio is Geographic.lib instead of
    GeographicLib.lib (compatible with makefiles).
  - Make classes NaN aware.
  - Use cell arrays for MGRS strings in Matlab.
  - Add solution/project files for Visual Studio 2010 (32-bit only).
  - Use C++11 static_assert and math functions, if available.

 - <a href="http://geographiclib.sf.net/1.6">Version 1.6</a>
   (released 2010-11-23)
  - FIX BUG introduced in GeographicLib::Geoid in version 1.5 (found by
    Dave Edwards).

 - <a href="http://geographiclib.sf.net/1.5">Version 1.5</a>
   (released 2010-11-19)
  - Improve area calculations for small polygons.
  - Add -s and -r flags to <a href="Planimeter.1.html">Planimeter</a>.
  - Improve the accuracy of GeographicLib::LambertConformalConic using
    divided differences.
  - FIX BUG in meridian convergence returned by
    LambertConformalConic::Forward.
  - Add optional threadsafe parameter to GeographicLib::Geoid
    constructor.  WARNING: This changes may break binary compatibility
    with previous versions of %GeographicLib.  However, the library is
    source compatible.
  - Add GeographicLib::OSGB.
  - Matlab and Octave interfaces to GeographicLib::UTMUPS,
    GeographicLib::MGRS, GeographicLib::Geoid, GeographicLib::Geodesic
    provided.
  - Minor changes
    - explicitly turn on optimization in Visual Studio 2008 projects;
    - add missing dependencies in some Makefiles;
    - move pi() and degree() from GeographicLib::Constants to
      GeographicLib::Math;
    - introduce GeographicLib::Math::extended type to aid testing;
    - add GeographicLib::Math::epi() and GeographicLib::Math::edegree().
    - fixes to compile under cygwin;
    - tweak expression used to find latitude from conformal latitude.

 - <a href="http://geographiclib.sf.net/1.4">Version 1.4</a>
   (released 2010-09-12)
  - Changes to GeographicLib::Geodesic and GeographicLib::GeodesicLine:
    - FIX BUG in Geodesic::Inverse with prolate ellipsoids;
    - add area computations to Geodesic::Direct and Geodesic::Inverse;
    - add geodesic areas to geodesic test set;
    - make GeodesicLine constructor public;
    - change longitude series in Geodesic into Helmert-like form;
    - ensure that equatorial geodesics have cos(alpha0) = 0 identically;
    - generalize interface for Geodesic and GeodesicLine;
    - split GeodesicLine and Geodesic into different files;
    - signal convergence failure in Geodesic::Inverse with NaNs;
    - deprecate one function in Geodesic and two functions in
      GeodesicLine;
    - deprecate -n option for <a href="GeodSolve.1.html">Geod</a>.
    .
    WARNING: These changes may break binary compatibility with previous
    versions of %GeographicLib.  However, the library is source
    compatible (with the proviso that GeographicLib/GeodesicLine.hpp may
    now need to be included).
  - Add the <a href="Planimeter.1.html">Planimeter</a> utility for
    computing the areas of geodesic polygons.
  - Improve iterative solution of GeographicLib::Gnomonic::Reverse.
  - Add GeographicLib::Geoid::ConvertHeight.
  - Add -msltohae, -haetomsl, and -z options to
    <a href="GeoidEval.1.html">GeoidEval</a>.
  - Constructors check that minor radius is positive.
  - Add overloaded Forward and Reverse functions to the projection
    classes which don't return the convergence (or azimuth) and scale.
  - Document function parameters and return values consistently.

 - <a href="http://geographiclib.sf.net/1.3">Version 1.3</a>
   (released 2010-07-21)
  - Add GeographicLib::Gnomonic, the ellipsoid generalization of the
    gnomonic projection.
  - Add -g and -e options to
    <a href="GeodesicProj.1.html">EquidistantTest</a>.
  - Use fixed-point notation for output from
    <a href="CartConvert.1.html">CartConvert</a>,
    <a href="GeodesicProj.1.html">EquidistantTest</a>,
    <a href="TransverseMercatorProj.1.html">TransverseMercatorTest</a>.
  - PolarStereographic:
    - Improved conversion to conformal coordinates;
    - Fix bug with scale at opposite pole;
    - Complain if latitude out of range in SetScale.
  - Add GeographicLib::Math::NaN().
  - Add long double version of hypot for Windows.
  - Add EllipticFunction::E(real).
  - Update references to Geotrans in MGRS documentation.
  - Speed up tmseries.mac.

 - <a href="http://geographiclib.sf.net/1.2">Version 1.2</a>
   (released 2010-05-21)
  - FIX BUGS in GeographicLib::Geodesic,
    - wrong azimuth returned by Direct if point 2 is on a pole;
    - Inverse sometimes fails with very close points.
  - Improve calculation of scale in GeographicLib::CassiniSoldner,
    - add GeographicLib::GeodesicLine::Scale,
      GeographicLib::GeodesicLine::EquatorialAzimuth, and
      GeographicLib::GeodesicLine::EquatorialArc;
    - break friend connection between CassiniSoldner and Geodesic.
  - Add GeographicLib::DMS::DecodeAngle and
    GeographicLib::DMS::DecodeAzimuth.  Extend
    GeographicLib::DMS::Decode and GeographicLib::DMS::Encode to deal
    with distances.
  - Code and documentation changes in GeographicLib::Geodesic and
    GeographicLib::Geocentric for consistency with
    the forthcoming paper on geodesics.
  - Increase order of series using in GeographicLib::Geodesic to 6 (full
    accuracy maintained for ellipsoid flattening < 0.01).
  - Macro __NO_LONG_DOUBLE_MATH to disable use of long double.
  - Correct declaration of GeographicLib::Math::isfinite to return a bool.
  - Changes in the \ref utilities,
    - improve error reporting when parsing command line arguments;
    - accept latitudes and longitudes in decimal degrees or degrees,
      minutes, and seconds, with optional hemisphere designators;
    - <a href="GeoConvert.1.html">GeoConvert</a> -z accepts zone or
      zone+hemisphere;
    - <a href="GeoidEval.1.html">GeoidEval</a> accepts any of the input
      formats used by <a href="GeoConvert.1.html">GeoConvert</a>;
    - <a href="CartConvert.1.html">CartConvert</a> allows the ellipsoid
      to be specified with -e.

 - <a href="http://geographiclib.sf.net/1.1">Version 1.1</a>
   (released 2010-02-09)
  - FIX BUG (introduced in 2009-03) in EllipticFunction::E(sn,cn,dn).
  - Increase accuracy of scale calculation in TransverseMercator and
    TransverseMercatorExact.
  - Code and documentation changes for consistency with
    <a  href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>

 - <a href="http://geographiclib.sf.net/1.0">Version 1.0</a>
   (released 2010-01-07)
  - Add autoconf configuration files.
  - BUG FIX: Improve initial guess for Newton's method in
    PolarStereographic::Reverse.  (Previously this failed to converge
    when the co-latitude exceeded about 130 deg.)
  - Constructors for TransverseMercator, TransverseMercatorExact,
    PolarStereographic, Geocentric, and Geodesic now check for obvious
    problems with their arguments and throw an exception if necessary.
  - Most classes now include inspector functions such as MajorRadius()
    so that you can determine how instances were constructed.
  - Add GeographicLib::LambertConformalConic class.
  - Add GeographicLib::PolarStereographic::SetScale to allow the
    latitude of true scale to be specified.
  - Add solution and project files for Visual Studio 2008.
  - Add GeographicLib::GeographicErr for exceptions.
  - GeographicLib::Geoid changes:
    - BUG FIX: fix typo in GeographicLib::Geoid::Cache which could cause
      a segmentation fault in some cases when the cached area spanned
      the prime meridian.
    - Include sufficient edge data to allow heights to be returned for
      cached area without disk reads;
    - Add inspector functions to query the extent of the cache.

 - <a href="http://geographiclib.sf.net/2009-11">Version 2009-11</a>
   (released 2009-11-03)
  - Allow specification of "closest UTM zone" in GeographicLib::UTMUPS
    and <a href="GeoConvert.1.html">GeoConvert</a> (via -t option).
  - Utilities now complain is there are too many tokens on input lines.
  - Include real-to-real versions of GeographicLib::DMS::Decode and
    GeographicLib::DMS::Encode.
  - More house-cleaning changes:
    - Ensure that functions which return results through reference
      arguments do not alter the arguments when an exception is thrown.
    - Improve accuracy of GeographicLib::MGRS::Forward.
    - Include more information in some error messages.
    - Improve accuracy of inverse hyperbolic functions.
    - Fix the way GeographicLib::Math functions handle different precisions.

 - <a href="http://geographiclib.sf.net/2009-10">Version 2009-10</a>
   (released 2009-10-18)
  - Change web site to http://geographiclib.sourceforge.net
  - Several house-cleaning changes:
    - Change from the a flat directory structure to a more easily
      maintained one.
    - Introduce Math class for common mathematical functions (in
      Constants.hpp).
    - Use Math::real as the type for all real quantities.  By default this
      is typedef'ed to double; and the library should be installed this
      way.
    - Eliminate const reference members of AzimuthalEquidistant,
      CassiniSoldner and LocalCartesian so that they may be copied.
    - Make several constructors explicit.  Disallow some constructors.
      Disallow copy constructor/assignment for Geoid.
    - Document least squares formulas in Geoid.cpp.
    - Use unsigned long long for files positions of geoid files in Geoid.
    - Introduce optional mgrslimits argument in UTMUPS::Forward and
      UTMUPS::Reverse to enforce stricter MGRS limits on eastings and
      northings.
    - Add 64-bit targets in Visual Studio project files.

 - <a href="http://geographiclib.sf.net/2009-09">Version 2009-09</a>
   (released 2009-09-01)
  - Add GeographicLib::Geoid and
    <a href="GeoidEval.1.html">GeoidEval</a> utility.

 - <a href="http://geographiclib.sf.net/2009-08">Version 2009-08</a>
   (released 2009-08-14)
  - Add GeographicLib::CassiniSoldner class and
    <a href="GeodesicProj.1.html">EquidistantTest</a> utility.
  - Fix bug in GeographicLib::Geodesic::Inverse where NaNs were
    sometimes returned.
  - INCOMPATIBLE CHANGE: AzimuthalEquidistant now returns the reciprocal
    of the azimuthal scale instead of the reduced length.
  - Add -n option to <a href="GeoConvert.1.html">GeoConvert</a>.

 - <a href="http://geographiclib.sf.net/2009-07">Version 2009-07</a>
   (released 2009-07-16)
  - Speed up the series inversion code in tmseries.mac and geod.mac.
  - Reference Borkowski in section on \ref geocentric.

 - <a href="http://geographiclib.sf.net/2009-06">Version 2009-06</a>
   (released 2009-06-01)
  - Add routines to decode and encode zone+hemisphere to GeographicLib::UTMUPS.
  - Clean up code in GeographicLib::Geodesic.

 - <a href="http://geographiclib.sf.net/2009-05">Version 2009-05</a>
   (released 2009-05-01)
  - Improvements to GeographicLib::Geodesic:
    - more economical series expansions,
    - return reduced length (as does the
      <a href="GeodSolve.1.html">Geod</a> utility),
    - improved calculation of starting point for inverse method,
    - use reduced length to give derivative for Newton's method.
  - Add GeographicLib::AzimuthalEquidistant class.
  - Make GeographicLib::Geocentric, GeographicLib::TransverseMercator,
    and GeographicLib::PolarStereographic classes work with prolate
    ellipsoids.
  - <a href="CartConvert.1.html">CartConvert</a> checks its inputs more
    carefully.
  - Remove reference to defunct Constants.cpp from GeographicLib.vcproj.

 - <a href="http://geographiclib.sf.net/2009-04">Version 2009-04</a>
   (released 2009-04-01)
  - Use compile-time constants to select the order of series in
    GeographicLib::TransverseMercator.
  - 2x unroll of Clenshaw summation to avoid data shuffling.
  - Simplification of GeographicLib::EllipticFunction::E.
  - Use STATIC_ASSERT for compile-time checking of constants.
  - Improvements to GeographicLib::Geodesic:
    - compile-time option to change order of series used,
    - post Maxima code for generating the series,
    - tune the order of series for double,
    - improvements in the selection of starting points for Newton's
      method,
    - accept and return spherical arc lengths,
    - works with both oblate and prolate ellipsoids,
    - add -a, -e, -b options to the <a href="GeodSolve.1.html">Geod</a>
      utility.

 - <a href="http://geographiclib.sf.net/2009-03">Version 2009-03</a>
   (released 2009-03-01)
  - Add GeographicLib::Geodesic and the
    <a href="GeodSolve.1.html">Geod</a> utility.
  - Declare when no exceptions are thrown by functions.
  - Minor changes to GeographicLib::DMS class.
  - Use invf = 0 to mean a sphere in constructors to some classes.
  - The makefile creates a library and includes an install target.
  - Rename GeographicLib::ECEF to GeographicLib::Geocentric, ECEFConvert
    to <a href="CartConvert.1.html">CartConvert</a>.
  - Use inline functions to define constant doubles in Constants.hpp.

 - <a href="http://geographiclib.sf.net/2009-02">Version 2009-02</a>
   (released 2009-01-30)
  - Fix documentation of constructors (flattening -> inverse
    flattening).
  - Use std versions of math functions.
  - Add GeographicLib::ECEF and GeographicLib::LocalCartesian classes
    and the ECEFConvert utility.
  - Gather the documentation on the \ref utilities onto one page.

 - <a href="http://geographiclib.sf.net/2009-01">Version 2009-01</a>
   (released 2009-01-12)
  - First proper release of library.
  - More robust GeographicLib::TransverseMercatorExact:
    - Introduce \e extendp version of constructor,
    - Test against extended test data,
    - Optimize starting positions for Newton's method,
    - Fix behavior near all singularities,
    - Fix order dependence in C++ start-up code,
    - Improved method of computing scale and convergence.
  - Documentation on transverse Mercator projection.
  - Add GeographicLib::MGRS, GeographicLib::UTMUPS, etc.

 - Version 2008-09
  - Ad hoc posting of information on the transverse Mercator projection.

<center>
Back to \ref geocentric.  Up to \ref contents.
</center>
**********************************************************************/
