12345qiupeng
/
gtsam
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
gtsam/gtsam/3rdparty/GeographicLib/doc/GeographicLib.dox

5546 lines
248 KiB
Plaintext
Raw Blame History

// -*- text -*-
/**
* \file GeographicLib.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 GeographicLib library
\author Charles F. F. Karney (charles@karney.com)
\version 1.34
\date 2013-12-11
\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. Also included is a .NET
wrapper library <a href="NET/index.html">NETGeographicLib</a>
which exposes the functionality to .NET applications. 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.34.tar.gz">
GeographicLib-1.34.tar.gz</a>
- <a href="http://sf.net/projects/geographiclib/files/distrib/GeographicLib-1.34.zip">
GeographicLib-1.34.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 -b r1.34 git://git.code.sf.net/p/geographiclib/code geographiclib
\endverbatim
There are also binary installers for Windows available at
- <a href="http://sf.net/projects/geographiclib/files/distrib/GeographicLib-1.34-win32.exe">
GeographicLib-1.34-win32.exe</a>
- <a href="http://sf.net/projects/geographiclib/files/distrib/GeographicLib-1.34-win64.exe">
GeographicLib-1.34-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), 11
(2012), and 12 (2013) compiled for 32 bit and 64 bit.
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.34.tar.gz">
GeographicLib-1.34.tar.gz</a> or
<a href="http://sourceforge.net/projects/geographiclib/files/distrib/GeographicLib-1.34.zip">
GeographicLib-1.34.zip</a> (or
<a href="http://sourceforge.net/projects/geographiclib/files/distrib/GeographicLib-1.34-win32.exe">
GeographicLib-1.34-win32.exe</a> or
<a href="http://sourceforge.net/projects/geographiclib/files/distrib/GeographicLib-1.34-win64.exe">
GeographicLib-1.34-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 library for use by Visual Studio 10.)
- \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). The minimum version of cmake supported for building
%GeographicLib is 2.8.4 (which was released on 2011-02-16). (Actually,
a few earlier versions will also work for non-Windows platforms.)
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.34.tar.gz
unzip -q GeographicLib-1.34.zip \endverbatim
then enter the directory created with one of \verbatim
cd GeographicLib-1.34 \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.34 ..
cmake -G "Visual Studio 9 2008" -D CMAKE_INSTALL_PREFIX=C:/pkg-vc9/GeographicLib-1.34 ..
\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.34. 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;
- python interface: OFF: python; ON: lib/python/site-packages;
- matlab interface: OFF: matlab; ON: libexec/GeographicLib/matlab;
.
- <code>CMAKE_INSTALL_PREFIX</code> (default: <code>/usr/local</code>
on non-Windows systems, <code>C:/Program Files/GeographicLib</code>
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.34
..</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>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>GEOGRAPHICLIB_LIB_TYPE</code> (allowed values: SHARED, STATIC, or
BOTH), specifies the types of libraries build. The default is
STATIC for Windows and SHARED otherwise. If building %GeographicLib
for sytem-wide use, BOTH is recommended, because this provides users
with the choice of which library to use.
- <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_DOCUMENTATION</code> (default: OFF). 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>BUILD_NETGEOGRAPHICLIB</code> (default: OFF). If set to ON,
build the managed C++ wrapper library
<a href="NET/index.html">NETGeographicLib</a>. This only makes
sense for Windows systems.
- <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
Possible additional targets are \verbatim
make matlabinterface (only for the Release configuration on Windows)
make dist
make exampleprograms
make netexamples (supported only for Release configuration) \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
For maximum flexibility, it's a good idea to build and install both
the Debug and Release versions of the library (in that order). The
installation directories will then contain the release versions of the
tools and <i>both</i> versions (debug and release) of the libraries.
If you use cmake to configure and build your programs, then the right
version of the library (debug vs. release) will automatically be used.
- 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.34.tar.gz \endverbatim
then enter the directory created \verbatim
cd GeographicLib-1.34 \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.34.tar.gz \endverbatim
then enter the directory created \verbatim
cd GeographicLib-1.34 \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.34.zip \endverbatim
- Open GeographicLib-1.34/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/GeographicLib.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.
The Visual Studio 10 solution also contains projects to build
<a href="NET/index.html">NETGeographicLib</a> and the C# example project.
\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 and shared versions of library are also
provided; these can only be used by Visual Studio 2010 (in either
release or debug mode). However, if you plan to use the library, it may
be 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.34-win32.exe">
GeographicLib-1.34-win32.exe</a> or
<a href="http://sourceforge.net/projects/geographiclib/files/distrib/GeographicLib-1.34-win64.exe">
GeographicLib-1.34-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.34 or C:\\pkg-vc10-x64\\GeographicLib-1.34,
- 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.34/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 libraries were built using
using Visual Studio 10 in release and debug modes. The utilities were
linked against the release-mode shared library. The Matlab interface
was compiled with Matlab R2013a 64-bit, however it may work with some
other 64-bit versions of Matlab.
<a href="NET/index.html">NETGeographicLib</a> library dlls (release and
debug) are included with the binary installers; these are linked against
the shared library for %GeographicLib.
\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 GeographicLib.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.34.
- <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.34.
.
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.34.tar.gz
GeographicLib-1.34.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>. On Windows, you should
configure %GeographicLib with <code>PACKAGE_DEBUG_LIBS</code>=ON, build both
Release and Debug versions of the library and finally build PACKAGE in
Release mode. This will get the release and debug versions of the
library included in the package. For example, the 64-bit binary
installer is created with \verbatim
cmake -G "Visual Studio 10 Win64" \
-D GEOGRAPHICLIB_LIB_TYPE=BOTH \
-D PACKAGE_DEBUG_LIBS=ON \
-D MATLAB_COMPILER=mex \
-D BUILD_NETGEOGRAPHICLIB=ON \
..
cmake --build . --config Debug --target ALL_BUILD
cmake --build . --config Release --target ALL_BUILD
cmake --build . --config Release --target matlabinterface
cmake --build . --config Release --target PACKAGE \endverbatim
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.34.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.
- Finally compile and link your code. You have two options here.
- Use cmake to build your package. If you are familiar with cmake
this typically will be far the simplest option.
- Set the include paths and linking options "manually".
- Building your code with cmake. In brief, the necessary steps are:
- include in your CMakeLists.txt files \verbatim
find_package (GeographicLib 1.34 REQUIRED)
include_directories (${GeographicLib_INCLUDE_DIRS})
add_definitions (${GeographicLib_DEFINITIONS})
add_executable (program source1.cpp source2.cpp)
target_link_libraries (program ${GeographicLib_LIBRARIES}) \endverbatim
(The <code>add_definitions</code> line is only needed for Windows
and can be omitted if you're using cmake version 2.8.11 or later.)
- configure your package, e.g., with \verbatim
mkdir BUILD
cd BUILD
cmake -G "Visual Studio 10" \
-D CMAKE_PREFIX_PATH=C:/pkg-vc10 \
-D CMAKE_PREFIX_PATH=C:/pkg-vc10/testgeographic \
.. \endverbatim
Note that you almost always want to configure and build your code
somewhere other than the source directory (in this case, we use the
BUILD subdirectory).
- build your package. On Linux and MacOS this usually involves just
running make. On Windows, you can load the solution file created by
cmake into Visual Studio; alternatively, you can get cmake to run
build your code with \verbatim
cmake --build . --config Release --target ALL_BUILD \endverbatim
You might also want to install your package (using "make install" or
build the "INSTALL" target with the command above).
.
The most import step is the find_package command. The cmake
documentation describes the locations searched by find_package (the
appropriate rule for %GeographicLib are those for "Config" mode lookups).
In brief, the locations that are searched are (from least specific to
most specific, i.e., in <i>reverse</i> order) are
- under the system paths, i.e., locations such as <code>C:/Program
Files</code> and <code>/usr/local</code>);
- frequently, it's necessary to search within a "package directory"
(or set of directories) for external dependencies; this is given by
a (semicolon separated) list of directories specified by the cmake
variable <code>CMAKE_PREFIX_PATH</code> (illustrated above);
- the package directory for %GeographicLib can be overridden with the
<i>environment variable</i> <code>GeographicLib_DIR</code> (which is the
directory under which %GeographicLib is installed);
- finally, if you need to point to a particular build of %GeographicLib,
define the <i>cmake variable</i> <code>GeographicLib_DIR</code>, which
specifies the directory containing the configuration file
<code>geographiclib-config.cmake</code> (for debugging this be the
top-level <i>build</i> directory, as opposed to <i>installation</i>
directory, for %GeographicLib).
.
Typically, specifying nothing or <code>CMAKE_PREFIX_PATH</code>
suffices. However the two <code>GeographicLib_DIR</code> variables allow
for a specific version to be chosen. On Windows systems (with Visual
Studio), find_package will only find 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.)
If %GeographicLib is found, then the following cmake variables are set:
- <code>GeographicLib_FOUND</code> = 1
- <code>GeographicLib_VERSION</code> = 1.34
- <code>GeographicLib_INCLUDE_DIRS</code>
- <code>GeographicLib_LIBRARIES</code> = one of the following two:
- <code>GeographicLib_SHARED_LIBRARIES</code> = %GeographicLib
- <code>GeographicLib_STATIC_LIBRARIES</code> = GeographicLib_STATIC
- <code>GeographicLib_DEFINITIONS</code> = one of the following two:
- <code>GeographicLib_SHARED_DEFINITIONS</code> = -DGEOGRAPHICLIB_SHARED_LIB=1
- <code>GeographicLib_STATIC_DEFINITIONS</code> = -DGEOGRAPHICLIB_SHARED_LIB=0
- <code>GeographicLib_LIBRARY_DIRS</code>
- <code>GeographicLib_BINARY_DIRS</code>
.
Either of <code>GeographicLib_SHARED_LIBRARIES</code> or
<code>GeographicLib_STATIC_LIBRARIES</code> may be empty, if that version
of the library is unavailable. If you require a specific version,
SHARED or STATIC, of the library, add a <code>COMPONENTS</code> clause
to find_package, e.g.,
\verbatim
find_package (GeographicLib 1.34 REQUIRED COMPONENTS SHARED) \endverbatim
causes only packages which include the shared library to be found. If
the package includes both versions of the library, then
<code>GeographicLib_LIBRARIES</code> and
<code>GeographicLib_DEFINITIONS</code> are set to the shared versions,
unless you include \verbatim
set (GeographicLib_USE_STATIC_LIBS ON) \endverbatim
<i>before</i> the find_package command. You can check whether
<code>GeographicLib_LIBRARIES</code> refers to the shared or static
library with \verbatim
get_target_property(_LIBTYPE ${GeographicLib_LIBRARIES} TYPE) \endverbatim
which results in <code>_LIBTYPE</code> being set to
<code>SHARED_LIBRARY</code> or <code>STATIC_LIBRARY</code>.
On Windows, cmake takes care of linking to the release or debug
version of the library as appropriate. (This assumes that the Release
and Debug versions of the libraries were built and installed. This is
true for the Windows binary installer for %GeographicLib version 1.34 and
later.)
- Here are the steps to compile and link your code using %GeographicLib
"manually".
- 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 -funroll-loops -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
- If using the shared (or static) library with Visual Studio, define
the macro <code>GEOGRAPHICLIB_SHARED_LIB=1</code> (or
<code>0</code>), e.g.,
\verbatim
C/C++ -> Preprocessor -> Preprocessor Definitions = GEOGRAPHICLIB_SHARED_LIB=1
\endverbatim
This is only needed for Windows systems when both shared and static
libraries have been installed. (If you configure your package with
cmake, this definition is added automatically.)
- 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-i.lib (for shared library)
Linker -> Input -> Additional Dependencies = Geographic.lib (for static library)
Linker -> General -> Additional Library Directories = C:\pkg-vc10\Geographic\lib
\endverbatim
Note that the library name is <b>Geographic</b> and not
%GeographicLib. For the Debug version of your program on Windows
add "_d" to the library, e.g., Geographic_d-i.lib or
Geographic_d.lib.
- Tell the runtime environment where to find the shared library
(assuming you compiled %Geographic 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 or Geographic_d.dll is in the same directory as
your executable or else include the directory containing the dll in
your <code>PATH</code>.
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
Here is a complete CMakeList.txt files you can use to build this test
code using the installed library: \verbatim
project (geodesictest)
cmake_minimum_required (VERSION 2.8.4)
find_package (GeographicLib 1.34 REQUIRED)
if (NOT MSVC)
set (CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
endif ()
include_directories (${GeographicLib_INCLUDE_DIRS})
add_definitions (${GeographicLib_DEFINITIONS})
add_executable (${PROJECT_NAME} example-Geodesic-small.cpp)
target_link_libraries (${PROJECT_NAME} ${GeographicLib_LIBRARIES})
if (MSVC)
get_target_property (_LIBTYPE ${GeographicLib_LIBRARIES} TYPE)
if (_LIBTYPE STREQUAL "SHARED_LIBRARY")
# On Windows systems, copy the shared library to build directory
add_custom_command (TARGET ${PROJECT_NAME} POST_BUILD
COMMAND ${CMAKE_COMMAND} -E
copy $<TARGET_FILE:${GeographicLib_LIBRARIES}> ${CMAKE_CFG_INTDIR}
COMMENT "Copying shared library for GeographicLib")
endif ()
endif () \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).
An attempt is made to maintain backward compatibility with
%GeographicLib (especially at the level of your source code). Sometimes
it's necessary to take actions that depend on what version of
%GeographicLib is being used; for example, you may want to use a new
feature of %GeographicLib, but want your code still to work with older
versions. In that case, you can test the values of the macros
GEOGRAPHICLIB_VERSION_MAJOR, GEOGRAPHICLIB_VERSION_MINOR, and
GEOGRAPHICLIB_VERSION_PATCH; these expand to numbers (and the last one
is usually 0); these macros appeared starting in version 1.34. There's
also a macro GEOGRAPHICLIB_VERSION_STRING which expands to, e.g.,
"1.34"; this macro has been defined since version 1.9.
<center>
Back to \ref utilities. Forward to \ref other. Up to \ref contents.
</center>
**********************************************************************/
/**
\page other Implementations in other languages
<center>
Back to \ref organization. Forward to \ref geoid. Up to \ref contents.
</center>
Implementations of subsets of %GeographicLib are available in other languages
- \ref c-fortran
- \ref java.
- \ref javascript.
- \ref python.
- \ref matlab.
- \ref maxima.
- \ref dotnet.
\section c-fortran C and Fortran implementation
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 java Java implementation
The directory <code>java</code> contains implementations of
GeographicLib::Geodesic, GeographicLib::GeodesicLine, and
GeographicLib::PolygonArea in Java. This implementation is entirely
self-contained and does 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="java/index.html">Java library for geodesics</a>.
\section javascript JavaScript implementation
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 implementation
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
Another illustrative exercise is finding the point midway between JFK
Airport to Singapore Changi Airport
\code
from geographiclib.geodesic import Geodesic
# Coordinates of airports
lat1, lon1 = 40.640, -73.779 # JFK
lat2, lon2 = 1.359, 103.989 # SIN
# Compute path from 1 to 2
g = Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2)
# Compute midpoint starting at 1
h = Geodesic.WGS84.Direct(lat1, lon1, g['azi1'], g['s12']/2)
print(h['lat2'], h['lon2']);
\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 implementations
The <code>matlab</code> 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.34/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 MEX-files:
Would you like mex to locate installed compilers [y]/n? y
Select a compiler:
[1] Microsoft Visual C++ 2012 in C:\Program Files (x86)\Microsoft Visual Studio 11.0
[2] Microsoft Visual C++ 2010 in C:\Program Files (x86)\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 ..
cmake --build . --config Release --target matlabinterface \endverbatim
(Note that only the Release configuration is supported for Matlab.)
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 R2013a 64-bit).<br>
On Linux systems, you can compile the interface for use with octave
instead by using \verbatim
cmake -D MATLAB_COMPILER=mkoctfile ..
make matlabinterface \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.34/matlab'
help geographiclibinterface
geographiclibinterface('C:/pkg-vc10/GeographicLib-1.34');
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.34/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.
\section dotnet .NET wrapper
This is a comprehensive wrapper library, written and maintained by Scott
Heiman, which exposes all of the functionality of %GeographicLib to the
.NET family of languages. For documentation, see
- <a href="NET/index.html">NETGeographicLib .NET wrapper library</a>.
<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 uses 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 GEOGRAPHICLIB_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 uses 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 GEOGRAPHICLIB_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 uses 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
- \ref geodshort
.
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
- The wikipedia page,
<a href="http://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid">
Geodesics on an ellipsoid</a>.
\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.
\image html meridian-measures.png "Comparison of meridian measures"
\section geodshort Short geodesics
Here we describe Bowring's method for solving the inverse geodesic
problem in the limit of short geodesics and contrast it with the great
circle solution using Bessel's auxiliary sphere. References:
- B. R. Bowring, The Direct and Inverse Problems for Short Geodesic
Lines on the Ellipsoid, Surveying and Mapping 41(2), 135--141 (1981).
- R. H. Rapp,
<a href="http://hdl.handle.net/1811/24333">
Geometric Geodesy, Part I</a>, Ohio State Univ. (1991), Sec. 6.5.
Bowring considers the conformal mapping of the ellipsoid to a sphere of
radius \f$ R \f$ such that circles of latitude and meridians are
preserved (and hence the azimuth of a line is preserved). Let \f$
(\phi, \lambda) \f$ and \f$ (\phi', \lambda') \f$ be the latitude and
longitude on the ellipsoid and sphere respectively. Define isometric
latitudes for the sphere and the ellipsoid as
\f[
\begin{aligned}
\psi' &= \sinh^{-1} \tan \phi', \\
\psi &= \sinh^{-1} \tan \phi - e \tanh^{-1}(e \sin\phi).
\end{aligned}
\f]
The most general conformal mapping satisfying Bowring's conditions is
\f[
\psi' = A \psi + K, \quad \lambda' = A \lambda,
\f]
where \f$ A \f$ and \f$ K \f$ are constants. (In fact a constant can be
added to the equation for \f$ \lambda' \f$, but this does affect the
analysis.) The scale of this mapping is
\f[
m(\phi) = \frac{AR}{\nu}\frac{\cos\phi'}{\cos\phi},
\f]
where \f$ \nu = a/\sqrt{1 - e^2\sin^2\phi} \f$ is the transverse radius
of curvature. (Note that in Bowring's Eq. (10), \f$ \phi \f$ should be
replaced by \f$ \phi' \f$.) The mapping from the ellipsoid to the sphere
depends on three parameters \f$ R, A, K \f$. These will be selected to
satisfy certain conditions at some representative latitude \f$ \phi_0
\f$. Two possible choices are given below.
\subsection bowring Bowring's method
Bowring (1981) requires that
\f[
m(\phi_0) = 1,\quad
\left.\frac{dm(\phi)}{d\phi}\right|_{\phi=\phi_0} = 0,\quad
\left.\frac{d^2m(\phi)}{d\phi^2}\right|_{\phi=\phi_0} = 0,
\f]
i.e, \f$m\approx 1\f$ in the vicinity of \f$\phi = \phi_0\f$.
This gives
\f[
\begin{aligned}
R &= \frac{\sqrt{1 + e'^2}}{B^2} a, \\
A &= \sqrt{1 + e'^2 \cos^4\phi_0}, \\
\tan\phi'_0 &= \frac1B \tan\phi_0,
\end{aligned}
\f]
where \f$ e' = e/\sqrt{1-e^2} \f$ is the second eccentricity, \f$ B =
\sqrt{1+e'^2\cos^2\phi_0} \f$, and \f$ K \f$ is defined implicitly by
the equation for \f$\phi'_0\f$. The radius \f$ R \f$ is the (Gaussian)
mean radius of curvature of the ellipsoid at \f$\phi_0\f$ (so near
\f$\phi_0\f$ the ellipsoid can be deformed to fit the sphere snugly).
The third derivative of \f$ m \f$ is given by
\f[
\left.\frac{d^3m(\phi)}{d\phi^3}\right|_{\phi=\phi_0} =
\frac{-2e'^2\sin2\phi_0}{B^4}.
\f]
The method for solving the inverse problem between two nearby points \f$
(\phi_1, \lambda_1) \f$ and \f$ (\phi_2, \lambda_2) \f$ is as follows:
Set \f$\phi_0 = (\phi_1 + \phi_2)/2\f$. Compute \f$ R, A, \phi'_0 \f$,
and hence find \f$ (\phi'_1, \lambda'_1) \f$ and \f$ (\phi'_2,
\lambda'_2) \f$. Finally, solve for the great circle on a sphere of
radius \f$ R \f$; the resulting distance and azimuths are good
approximations for the corresponding quantities for the ellipsoidal
geodesic.
Consistent with the accuracy of this method, we can compute
\f$\phi'_1\f$ and \f$\phi'_2\f$ using a Taylor expansion about
\f$\phi_0\f$. This also avoids numerical errors that arise from
subtracting nearly equal quantities when using the equation for
\f$\phi'\f$ directly. Write \f$\Delta \phi = \phi - \phi_0\f$ and
\f$\Delta \phi' = \phi' - \phi'_0\f$; then we have
\f[
\Delta\phi' \approx
\frac{\Delta\phi}B \biggl[1 +
\frac{\Delta\phi}{B^2}\frac{e'^2}2
\biggl(3\sin\phi_0\cos\phi_0 +
\frac{\Delta\phi}{B^2}
\bigl(B^2 - \sin^2\phi_0(2 - 3 e'^2 \cos^2\phi_0)\bigr)\biggr)\biggr],
\f]
where the error is \f$O(f\Delta\phi^4)\f$.
This is essentially Bowring's method. Significant differences between
this result, "Bowring (improved)", compared to Bowring's paper, "Bowring
(original)", are:
- Bowring elects to use \f$\phi_0 = \phi_1\f$. This simplifies the
calculations somewhat but increases the error by about a factor of
4.
- Bowring's expression for \f$ \Delta\phi' \f$ is only accurate in the
limit \f$ e' \rightarrow 0 \f$.
.
In fact, arguably, the highest order \f$O(f\Delta\phi^3)\f$ terms should
be dropped altogether. Their inclusion does result in a better estimate
for the distance. However, if your goal is to generate both accurate
distances \e and accurate azimuths, then \f$\Delta\phi\f$ needs to be
restricted sufficiently to allow these terms to be dropped to give the
"Bowring (truncated)" method.
With highly eccentric ellipsoids, the parametric latitude \f$ \beta \f$
is a better behaved independent variable to use. In this case, \f$
\phi_0 \f$ is naturally defined using \f$\beta_0 = (\beta_1 +
\beta_2)/2\f$ and in terms of \f$\Delta\beta = \beta - \beta_0\f$, we
have
\f[
\Delta\phi' \approx
\frac{\Delta\beta}{B'} \biggl[1 +
\frac{\Delta\beta}{B'^2}\frac{e'^2}2
\biggl(\sin\beta_0\cos\beta_0 +
\frac{\Delta\beta}{3B'^2}
\bigl( \cos^2\beta_0 - \sin^2\beta_0 B'^2\bigr)
\biggr)\biggr],
\f]
where \f$B' = \sqrt{1+e'^2\sin^2\beta_0} = \sqrt{1+e'^2}/B\f$, and the
error once again is \f$O(f\Delta\phi^4)\f$. This is the &quot;Bowring
(using \f$\beta\f$)&quot; method.
\subsection auxsphere Bessel's auxiliary sphere
%GeographicLib's uses the auxiliary sphere method of Legendre, Bessel,
and Helmert. For short geodesics, this is equivalent to picking
\f$ R, A, K \f$ so that
\f[
m(\phi_0) = 1,\quad
\left.\frac{dm(\phi)}{d\phi}\right|_{\phi=\phi_0} = 0,\quad
\tan\phi'_0 = (1 - f) \tan\phi_0.
\f]
Bowring's requirement that the second derivative of \f$m\f$ vanish has
been replaced by the last relation which states that \f$\phi'_0 =
\beta_0\f$, the parametric latitude corresponding to \f$\phi_0\f$. This
gives
\f[
\begin{aligned}
R &= B'(1-f)a, \\
A &= \frac1{B'(1-f)}, \\
\left.\frac{d^2m(\phi)}{d\phi^2}\right|_{\phi=\phi_0} &=
-e^2B'^2\sin^2\phi_0.
\end{aligned}
\f]
Similar to Bowring's method, we can compute \f$\phi'_1\f$ and
\f$\phi'_2\f$ using a Taylor expansion about \f$\beta_0\f$. This results
in the simple expression
\f[
\Delta\phi' \approx \Delta\beta,
\f]
where the error is \f$O(f\Delta\beta^2)\f$.
\subsection shorterr Estimating the accuracy
In assessing the accuracy of these methods we use two metrics:
- The absolute error in the distance.
- The consistency of the predicted azimuths. Imagine starting
ellipsoidal geodesics at \f$ (\phi_1, \lambda_1) \f$ and \f$ (\phi_2,
\lambda_2) \f$ with the predicted azimuths. What is the distance
between them when they are extended a distance \f$ a \f$ beyond the
second point?
.
(The second metric is much more stringent.) We may now compare the
methods by asking for a bound to the length of a geodesic which ensures
that the one or other of the errors fall below 1 mm (an "engineering"
definition of accurate) or 1 nm (1 nanometer, about the round-off
limit).
<center>
<table>
<caption>Maximum distance that can be used in various methods for
computing short geodesics while keeping the errors within prescribed
bounds</caption>
<tr>
<th rowspan="2">method
<th colspan="2"><center>distance metric</center></th>
<th colspan="2"><center>azimuth metric</center></th>
<tr>
<th>1 mm error</th>
<th>1 nm error</th>
<th>1 mm error</th>
<th>1 nm error</th>
<tr>
<td>Bowring (original)
<td><center>87 km</center>
<td><center>870 m</center>
<td><center>35 km</center>
<td><center>350 m</center>
<tr>
<td>Bowring (improved)
<td><center>180 km</center>
<td><center>1.8 km</center>
<td><center>58 km</center>
<td><center>580 m</center>
<tr>
<td>Bowring (truncated)
<td><center>52 km</center>
<td><center>520 m</center>
<td><center>52 km</center>
<td><center>520 m</center>
<tr>
<td>Bowring (using \f$\beta\f$)
<td><center>380 km</center>
<td><center>24 km</center>
<td><center>60 km</center>
<td><center>600 m</center>
<tr>
<td>Bessel's aux. sphere
<td><center>42 km</center>
<td><center>420 m</center>
<td><center>1.7 km</center>
<td><center>1.7 m</center>
</table>
</center>
For example, if you're only interested in measuring distances and an
accuracy of 1 mm is sufficient, then Bowring's improved method can be
used for distances up to 180 km. On the other hand, %GeographicLib uses
Bessel's auxiliary sphere and we require both the distance and the
azimuth to be accurate, so the great circle approximation can only be
used for distances less than 1.7 m. The reason that %GeographicLib does
not use Bowring's method is that the information necessary for auxiliary
sphere method is already available as part of the general solution and,
as much as possible, we allow all geodesics to be computed by the
general method.
<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
ellipsoid (with 3 unequal axes) can be reduced to 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.
See also
- The wikipedia page,
<a href="http://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid#Geodesics_on_a_triaxial_ellipsoid">
Geodesics on a triaxial ellipsoid</a>.
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="https://upload.wikimedia.org/wikipedia/commons/b/bd/Triaxial_ellipsoid_coordinate_system.svg"
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="https://upload.wikimedia.org/wikipedia/commons/a/a9/Circumpolar_geodesic_on_a_triaxial_ellipsoid_case_A.svg"
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="https://upload.wikimedia.org/wikipedia/commons/f/f3/Circumpolar_geodesic_on_a_triaxial_ellipsoid_case_B.svg"
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="https://upload.wikimedia.org/wikipedia/commons/8/83/Transpolar_geodesic_on_a_triaxial_ellipsoid_case_A.svg"
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="https://upload.wikimedia.org/wikipedia/commons/9/9c/Transpolar_geodesic_on_a_triaxial_ellipsoid_case_B.svg"
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 it
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="https://upload.wikimedia.org/wikipedia/commons/e/e5/Unstable_umbilical_geodesic_on_a_triaxial_ellipsoid.svg"
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>,
was 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&deg;, 90&deg;], lon
in [0&deg;, 90&deg;];
- 1000 entries randomly distributed on lat in [0&deg;, 90&deg;], lon =
0&deg;;
- 1000 entries randomly distributed on lat = 0&deg;, lon in [0&deg;,
90&deg;];
- 1000 entries randomly distributed on lat in [0&deg;, 90&deg;], lon =
90&deg;;
- 1000 entries close to lat = 90&deg; with lon in [0&deg;, 90&deg;];
- 1000 entries close to lat = 0&deg;, lon = 0&deg; with lat &ge;
0&deg;, lon &ge; 0&deg;;
- 1000 entries close to lat = 0&deg;, lon = 90&deg; with lat &ge;
0&deg;, lon &le; 90&deg;;
- 2000 entries close to lat = 0&deg;, lon = 90&deg; (1 &minus; \e e)
with lat &ge; 0&deg;;
- 25000 entries randomly distributed in lat in [&minus;89&deg;,
0&deg;], lon in [90&deg; (1 &minus; \e e), 90&deg;];
- 1000 entries randomly distributed on lat in [&minus;89&deg;, 0&deg;],
lon = 90&deg;;
- 1000 entries randomly distributed on lat in [&minus;89&deg;, 0&deg;],
lon = 90&deg; (1 &minus; \e e);
- 1000 entries close to lat = 0&deg;, lon = 90&deg; (lat < 0&deg;, lon
&le; 90&deg;);
- 1000 entries close to lat = 0&deg;, lon = 90&deg; (1 &minus; \e e)
(lat < 0&deg;, lon &le; 90&deg; (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.
\image html gauss-schreiber-graticule-a.png "Fig. 1(a)"
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.
\image html gauss-krueger-graticule-a.png "Fig. 1(b)"
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.
\image html thompson-tm-graticule-a.png "Fig. 1(c)"
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.
\image html gauss-krueger-error.png "Fig. 2"
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)
\image html gauss-krueger-graticule.png "Fig. 3(a)"
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;].
\image html gauss-krueger-convergence-scale.png "Fig. 3(b)"
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.
\image html thompson-tm-graticule.png "Fig. 3(c)"
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.34">Version 1.34</a>
(released 2013-12-11)
- Many changes in cmake support:
- minimum version of cmake needed increased to 2.8.4 (which was
released in 2011-02);
- allow building both shared and static librarys with <code>-D
GEOGRAPHICLIB_LIB_TYPE=BOTH</code>;
- both shared and static libraries (Release plus Debug) included in
binary installer;
- find_package uses COMPONENTS and GeographicLib_USE_STATIC_LIBS to
select the library to use;
- find_package version checking allows nmake and Visual Studio
generators to interoperate on Windows;
- find_package (%GeographicLib ...) requires that %GeographicLib be
capitalized correctly;
- on Unix/Linux, don't include the version number in directory for
the cmake configuration files;
- defaults for GEOGRAPHICLIB_DOCUMENTATION and
BUILD_NETGEOGRAPHICLIB are now OFF;
- the GEOGRAPHICLIB_EXAMPLES configuration parameter is no longer
used; cmake always configures to build the examples, but they are
not built by default (instead build targets: exampleprograms and
netexamples);
- matlab-all target renamed to matlabinterface;
- the configuration parameters PACKAGE_PATH and INSTALL_PATH are
now deprecated (use CMAKE_INSTALL_PREFIX instead);
- on Linux, the installed package is relocatable;
- on MacOSX, the installed utilities can find the shared library.
- Use a more precise value for GeographicLib::OSGB::CentralScale().
- Add Arc routines to python interface.
- The Geod utility has been removed; the same functionality lives on
with <a href="GeodSolve.1.html">GeodSolve</a> (introduced in
version 1.30).
- <a href="http://geographiclib.sf.net/1.33">Version 1.33</a>
(released 2013-10-08)
- Add <a href="NET/index.html">NETGeographic .NET wrapper library</a>
(courtesy of Scott Heiman).
- Make inspector functions in GeographicLib::Ellipsoid const.
- Add Accumulator.cpp to instantiate GeographicLib::Accumulator.
- Defer some of the initialization of GeographicLib::OSGB to when it
is first called.
- Fix bug in autoconf builds under MacOS.
- <a href="http://geographiclib.sf.net/1.32">Version 1.32</a>
(released 2013-07-12)
- Generalize C interface for polygon areas to allow vertices to be
specified incrementally.
- Fix way flags for C++11 support are determined.
- <a href="http://geographiclib.sf.net/1.31">Version 1.31</a>
(released 2013-07-01)
- Changes breaking binary compatibility (source compatibility is
maintained):
- overloaded versions of GeographicLib::DMS::Encode,
GeographicLib::EllipticFunction::EllipticFunction, and
GeographicLib::GeoCoords::DMSRepresentation, have been eliminated
by the use of optional arguments;
- correct the declaration of first arg to
GeographicLib::UTMUPS::DecodeEPSG.
- FIX BUG in GeographicLib::GravityCircle constructor (found by
Mathieu Peyr&eacute;ga) which caused bogus results for the gravity
disturbance and gravity anomaly vectors. (This only affected
calculations using GravityCircle. GravityModel calculations did
not suffer from this bug.)
- Improvements to the build:
- add macros GEOGRAPHICLIB_VERSION_{MAJOR,MINOR,PATCH} to Config.h;
- fix documentation for new version of perlpod;
- improving setting of runtime path for Unix-like systems with cmake;
- install PDB files when compiling with Visual Studio to aid
debugging;
- Windows binary release now uses Matlab R2013a (64-bit) and
uses the -largeArrayDims option.
- fixes to the way the Matlab interface routines are built (thanks
to Phil Miller and Chris F.).
- Changes to the geodesic routines:
- add \ref java of the geodesic routines (thanks to Skip Breidbach
for the maven support);
- FIX BUG: avoid altering input args in Fortran implementation;
- more systematic treatment of very short geodesic;
- fixes to python port so that they work with version 3.x, in
addition to 2.x (courtesy of Amato);
- accumulate the perimeter and area of polygons via a double-wide
accumulator in Fortran, C, and Matlab implementations (this is
already included in the other implementations);
- port GeographicLib::PolygonArea::AddEdge and
GeographicLib::PolygonArea::TestEdge to JavaScript and python
interfaces;
- include documentation on \ref geodshort.
- Unix scripts for downloading datasets,
geographiclib-get-{geoids,gravity,magnetic}, skip already download
models by default, unless the -f flag is given.
- FIX BUGS: meridian convergence and scale returned by
GeographicLib::TransverseMercatorExact was wrong at a pole.
- Improve efficiency of GeographicLib::MGRS::Forward by avoiding the
calculation of the latitude if possible (adapting an idea of Craig
Rollins).
- Fixes to the way the Matlab interface routines are built (thanks to
Phil Miller and Chris F.).
- <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 \ref triaxial 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>
**********************************************************************/
Reference in New Issue View Git Blame Copy Permalink