Major upgrade to geodesic support from Charles Karney (#197).

Syncs geodesic routines with GeographicLib.  Adds geodesic.3 man page. 
geod_* api exposed publically.  geodesic.h is installed.


git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@2333 4e78687f-474d-0410-85f9-8d5e500ac6b2

1 files changed, 493 insertions, 166 deletions

diff --git a/src/geodesic.h b/src/geodesic.h
index acf5d875..2c49d0dc 100644
--- a/src/geodesic.h
+++ b/src/geodesic.h
@@ -1,138 +1,123 @@
-/*
- * This is a C implementation of the geodesic algorithms described in
- *
- *   C. F. F. Karney,
- *   Algorithms for geodesics,
- *   J. Geodesy (2012);
- *   http://dx.doi.org/10.1007/s00190-012-0578-z
- *   Addenda: http://geographiclib.sf.net/geod-addenda.html
+/**
+ * \file geodesic.h
+ * \brief Header for the geodesic routines in C
  *
+ * This an implementation in C of the geodesic algorithms described in
+ * - C. F. F. Karney,
+ *   <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
+ *   Algorithms for geodesics</a>,
+ *   J. Geodesy <b>87</b>, 43--55 (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>.
+ * .
  * The principal advantages of these algorithms over previous ones (e.g.,
  * Vincenty, 1975) are
- *   * accurate to round off for abs(f) < 1/50;
- *   * the solution of the inverse problem is always found;
- *   * differential and integral properties of geodesics are computed.
+ * - accurate to round off for |<i>f</i>| &lt; 1/50;
+ * - the solution of the inverse problem is always found;
+ * - differential and integral properties of geodesics are computed.
  *
- * The shortest path between two points on the ellipsoid at (lat1, lon1) and
- * (lat2, lon2) is called the geodesic.  Its length is s12 and the geodesic
- * from point 1 to point 2 has forward azimuths azi1 and azi2 at the two end
- * points.
+ * The shortest path between two points on the ellipsoid at (\e lat1, \e
+ * lon1) and (\e lat2, \e lon2) is called the geodesic.  Its length is
+ * \e s12 and the geodesic from point 1 to point 2 has forward azimuths
+ * \e azi1 and \e azi2 at the two end points.
  *
  * Traditionally two geodesic problems are considered:
- *   * the direct problem -- given lat1, lon1, s12, and azi1, determine
- *     lat2, lon2, and azi2.
- *   * the inverse problem -- given lat1, lon1, lat2, lon2, determine
- *     s12, azi1, and azi2.
+ * - the direct problem -- given \e lat1, \e lon1, \e s12, and \e azi1,
+ *   determine \e lat2, \e lon2, and \e azi2.  This is solved by the function
+ *   geod_direct().
+ * - the inverse problem -- given \e lat1, \e lon1, \e lat2, \e lon2, determine
+ *   \e s12, \e azi1, and \e azi2.  This is solved by the function
+ *   geod_inverse().
  *
- * The ellipsoid is specified by its equatorial radius a (typically in meters)
- * and flattening f.  The routines are accurate to round off with double
- * precision arithmetic provided that abs(f) < 1/50; for the WGS84 ellipsoid,
- * the errors are less than 15 nanometers.  (Reasonably accurate results are
- * obtained for abs(f) < 1/5.)  Latitudes, longitudes, and azimuths are in
- * degrees.  Latitudes must lie in [-90,90] and longitudes and azimuths must
- * lie in [-540,540).  The returned values for longitude and azimuths are in
- * [-180,180).  The distance s12 is measured in meters (more precisely the same
- * units as a).
+ * The ellipsoid is specified by its equatorial radius \e a (typically in
+ * meters) and flattening \e f.  The routines are accurate to round off with
+ * double precision arithmetic provided that |<i>f</i>| &lt; 1/50; for the
+ * WGS84 ellipsoid, the errors are less than 15 nanometers.  (Reasonably
+ * accurate results are obtained for |<i>f</i>| &lt; 1/5.)  For a prolate
+ * ellipsoid, specify \e f &lt; 0.
  *
  * The routines also calculate several other quantities of interest
- *   * SS12 is the area between the geodesic from point 1 to point 2 and the
- *     equator; i.e., it is the area, measured counter-clockwise, of the
- *     quadrilateral with corners (lat1,lon1), (0,lon1), (0,lon2), and
- *     (lat2,lon2).  It is given in meters^2.
- *   * m12, the reduced length of the geodesic is defined such that if the
- *     initial azimuth is perturbed by dazi1 (radians) then the second point is
- *     displaced by m12 dazi1 in the direction perpendicular to the geodesic.
- *     m12 is given in meters.  On a curved surface the reduced length obeys a
- *     symmetry relation, m12 + m21 = 0.  On a flat surface, we have m12 = s12.
- *   * MM12 and MM21 are geodesic scales.  If two geodesics are parallel at
- *     point 1 and separated by a small distance dt, then they are separated by
- *     a distance MM12 dt at point 2.  MM21 is defined similarly (with the
- *     geodesics being parallel to one another at point 2).  MM12 and MM21 are
- *     dimensionless quantities.  On a flat surface, we have MM12 = MM21 = 1.
- *   * a12 is the arc length on the auxiliary sphere.  This is a construct for
- *     converting the problem to one in spherical trigonometry.  a12 is
- *     measured in degrees.  The spherical arc length from one equator crossing
- *     to the next is always 180 degrees.
- *
- * Simple interface
- *
- *    #include "geodesic.h"
- *
- *    double a, f, lat1, lon1, azi1, lat2, lon2, azi2, s12;
- *    struct Geodesic g;
- *
- *    GeodesicInit(&g, a, f);
- *    Direct(&g, lat1, lon1, azi1, s12, &lat2, &lon2, &azi2);
- *    Inverse(&g, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2);
- *
- * GeodesicInit initalizes g for the ellipsoid.  Subsequent calls to Direct and
- * Inverse solve the direct and inverse geodesic problems.
- *
- * Returning auxiliary quantities (continued from the previous example):
- *
- *    double a12, s12_a12, m12, M12, M21, S12;
- *    a12 = GenDirect(&g, lat1, lon1, azi1, 0, s12,
- *                    &lat1, &lat2, &azi2, 0, &m12, &M12, &M21, &S21);
- *    GenDirect(&g, lat1, lon1, azi1, 1, a12,
- *              &lat1, &lat2, &azi2, &s12, &m12, &M12, &M21, &S21);
- *    a12 = GenInverse(&g, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2,
- *                     &m12, &M12, &M21, &S12);
- *
- * GenDirect is a generalized version of Direct allowing the return of the
- * auxiliary quantities.  With the first variant (arcmode = 0), the length of
- * the geodesic is specified by the true length s12 and the arc length a12 is
- * returned as the function value.  In the second variant (arcmode = 1), the
- * length is specified by the arc length a12 (in degrees), and the true length
- * is returned via &s12.
- *
- *    a12 = GenInverse(&g, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2,
- *                     &m12, &M12, &M21, &S12);
- *
- * GenInverse is a generalized version of Inverse allowing the return of the
- * auxiliary quantities.
- *
- * Any of the "return" arguments &s12, etc. in these routines may be replaced
- * by 0, if you do not need some quantities computed.
- *
- * Computing multiple points on a geodesic.  This may be accomplished by
- * repeated invocations of Direct varying s12.  However, it is more efficient
- * to create a GeodesicLine object, as follows.
- *
- *    struct GeodesicLine l;
- *    int caps = 0;
+ * - \e S12 is the area between the geodesic from point 1 to point 2 and the
+ *   equator; i.e., it is the area, measured counter-clockwise, of the
+ *   quadrilateral with corners (\e lat1,\e lon1), (0,\e lon1), (0,\e lon2),
+ *   and (\e lat2,\e lon2).
+ * - \e m12, the reduced length of the geodesic is defined such that if
+ *   the initial azimuth is perturbed by \e dazi1 (radians) then the
+ *   second point is displaced by \e m12 \e dazi1 in the direction
+ *   perpendicular to the geodesic.  On a curved surface the reduced
+ *   length obeys a symmetry relation, \e m12 + \e m21 = 0.  On a flat
+ *   surface, we have \e m12 = \e s12.
+ * - \e M12 and \e M21 are geodesic scales.  If two geodesics are
+ *   parallel at point 1 and separated by a small distance \e dt, then
+ *   they are separated by a distance \e M12 \e dt at point 2.  \e M21
+ *   is defined similarly (with the geodesics being parallel to one
+ *   another at point 2).  On a flat surface, we have \e M12 = \e M21
+ *   = 1.
+ * - \e a12 is the arc length on the auxiliary sphere.  This is a
+ *   construct for converting the problem to one in spherical
+ *   trigonometry.  \e a12 is measured in degrees.  The spherical arc
+ *   length from one equator crossing to the next is always 180&deg;.
  *
- *    GeodesicLineInit(&l, &g, a, f, lat1, lon1, azi1, caps);
- *    Position(l, s12, &lat2, &lon2, &azi2)
+ * If points 1, 2, and 3 lie on a single geodesic, then the following
+ * addition rules hold:
+ * - \e s13 = \e s12 + \e s23
+ * - \e a13 = \e a12 + \e a23
+ * - \e S13 = \e S12 + \e S23
+ * - \e m13 = \e m12 \e M23 + \e m23 \e M21
+ * - \e M13 = \e M12 \e M23 &minus; (1 &minus; \e M12 \e M21) \e
+ *   m23 / \e m12
+ * - \e M31 = \e M32 \e M21 &minus; (1 &minus; \e M23 \e M32) \e
+ *   m12 / \e m23
  *
- * caps is a bit mask specifying the capabilities of the GeodesicLine object.
- * It should be an or'ed combination of
+ * The shortest distance returned by the solution of the inverse problem is
+ * (obviously) uniquely defined.  However, in a few special cases there are
+ * multiple azimuths which yield the same shortest distance.  Here is a
+ * catalog of those cases:
+ * - \e lat1 = &minus;\e lat2 (with neither at a pole).  If \e azi1 = \e
+ *   azi2, the geodesic is unique.  Otherwise there are two geodesics
+ *   and the second one is obtained by setting [\e azi1, \e azi2] = [\e
+ *   azi2, \e azi1], [\e M12, \e M21] = [\e M21, \e M12], \e S12 =
+ *   &minus;\e S12.  (This occurs when the longitude difference is near
+ *   &plusmn;180&deg; for oblate ellipsoids.)
+ * - \e lon2 = \e lon1 &plusmn; 180&deg; (with neither at a pole).  If
+ *   \e azi1 = 0&deg; or &plusmn;180&deg;, the geodesic is unique.
+ *   Otherwise there are two geodesics and the second one is obtained by
+ *   setting [\e azi1, \e azi2] = [&minus;\e azi1, &minus;\e azi2], \e
+ *   S12 = &minus;\e S12.  (This occurs when the \e lat2 is near
+ *   &minus;\e lat1 for prolate ellipsoids.)
+ * - Points 1 and 2 at opposite poles.  There are infinitely many
+ *   geodesics which can be generated by setting [\e azi1, \e azi2] =
+ *   [\e azi1, \e azi2] + [\e d, &minus;\e d], for arbitrary \e d.  (For
+ *   spheres, this prescription applies when points 1 and 2 are
+ *   antipodal.)
+ * - \e s12 = 0 (coincident points).  There are infinitely many geodesics
+ *   which can be generated by setting [\e azi1, \e azi2] = [\e azi1, \e
+ *   azi2] + [\e d, \e d], for arbitrary \e d.
  *
- *    LATITUDE        compute lat2 (in effect this is always set)
- *    LONGITUDE       compute lon2
- *    AZIMUTH         compute azi2 (in effect this is always set)
- *    DISTANCE        compute s12
- *    DISTANCE_IN     allow the length to be specified in terms of distance
- *    REDUCEDLENGTH   compute m12
- *    GEODESICSCALE   compute M12 and M21
- *    AREA            compute S12
- *    ALL             all of the above
+ * These routines are a simple transcription of the corresponding C++ classes
+ * in <a href="http://geographiclib.sf.net"> GeographicLib</a>.  The
+ * "class data" is represented by the structs geod_geodesic and
+ * geod_geodesicline and these are passed as initial arguments to the member
+ * functions.  (But note that the PolygonArea class has been replaced by a
+ * simple function interface and the area computation does use the accurate
+ * summing providing by the Accumulator class.)  Most of the internal comments
+ * have been retained.  However, in the process of transcription some
+ * documentation has been lost and the documentation for the C++ classes,
+ * GeographicLib::Geodesic, GeographicLib::GeodesicLine, and
+ * GeographicLib::PolygonArea, should be consulted.  The C++ code remains the
+ * "reference implementation".  Think twice about restructuring the internals
+ * of the C code since this may make porting fixes from the C++ code more
+ * difficult.
  *
- * caps = 0 is treated as LATITUDE | LONGITUDE | AZIMUTH | DISTANCE_IN (to
- * support the solution "standard" direct problem).
- *
- * There's also a generalized version of Position
- *
- *    a12 = GenPosition(&l, arcmode, s12_a12,
- *                      &lat2, &lon2, &azi2, &s12, &m12, &M12, &M21, &S12);
- *
- * See the documentation on GenDirect for the meaning of arcmode.
- *
- * Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2013) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * http://geographiclib.sourceforge.net/
  *
- * This file was distributed with GeographicLib 1.27.
- */
+ * This library was distributed with
+ * <a href="../index.html">GeographicLib</a> 1.30.
+ **********************************************************************/
 
 #if !defined(GEODESIC_H)
 #define GEODESIC_H 1
@@ -141,62 +126,404 @@
 extern "C" {
 #endif
 
-  struct Geodesic {
-    double a, f, f1, e2, ep2, n, b, c2, etol2;
+  /**
+   * The struct containing information about the ellipsoid.
+   **********************************************************************/
+  struct geod_geodesic {
+    double a;                   /**< the equatorial radius */
+    double f;                   /**< the flattening */
+    /**< @cond SKIP */
+    double f1, e2, ep2, n, b, c2, etol2;
     double A3x[6], C3x[15], C4x[21];
+    /**< @endcond */
   };
 
-  struct GeodesicLine {
-    double lat1, lon1, azi1;
-    double a, f, b, c2, f1, salp0, calp0, k2,
+  /**
+   * The struct containing information about a single geodesic.
+   **********************************************************************/
+  struct geod_geodesicline {
+    double lat1;                /**< the starting latitude */
+    double lon1;                /**< the starting longitude */
+    double azi1;                /**< the starting azimuth */
+    double a;                   /**< the equatorial radius */
+    double f;                   /**< the flattening */
+    /**< @cond SKIP */
+    double b, c2, f1, salp0, calp0, k2,
       salp1, calp1, ssig1, csig1, dn1, stau1, ctau1, somg1, comg1,
       A1m1, A2m1, A3c, B11, B21, B31, A4, B41;
     double C1a[6+1], C1pa[6+1], C2a[6+1], C3a[6], C4a[6];
-    unsigned caps;
+    /**< @endcond */
+    unsigned caps;              /**< the capabilities */
   };
 
-  void GeodesicInit(struct Geodesic* g, double a, double f);
-  void GeodesicLineInit(struct GeodesicLine* l,
-                        const struct Geodesic* g,
-                        double lat1, double lon1, double azi1, unsigned caps);
-
-  void Direct(const struct Geodesic* g,
-              double lat1, double lon1, double azi1, double s12,
-              double* plat2, double* plon2, double* pazi2);
-  void Inverse(const struct Geodesic* g,
-               double lat1, double lon1, double lat2, double lon2,
-               double* ps12, double* pazi1, double* pazi2);
-  void Position(const struct GeodesicLine* l, double s12,
-                double* plat2, double* plon2, double* pazi2);
-
-  double GenDirect(const struct Geodesic* g,
-                   double lat1, double lon1, double azi1,
-                   int arcmode, double s12_a12,
-                   double* plat2, double* plon2, double* pazi2,
-                   double* ps12, double* pm12, double* pM12, double* pM21,
-                   double* pS12);
-  double GenInverse(const struct Geodesic* g,
+  /**
+   * Initialize a geod_geodesic object
+   *
+   * @param[out] g the object to be initialized.
+   * @param[in] a the equatorial radius (meters).
+   * @param[in] f the flattening.
+   **********************************************************************/
+  void geod_init(struct geod_geodesic* g, double a, double f);
+
+  /**
+   * Initialize a geod_geodesicline object
+   *
+   * @param[out] l the object to be initialized.
+   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] lat1 latitude of point 1 (degrees).
+   * @param[in] lon1 longitude of point 1 (degrees).
+   * @param[in] azi1 azimuth at point 1 (degrees).
+   * @param[in] caps bitor'ed combination of geod_mask values specifying the
+   *   capabilities the geod_geodesicline object should possess, i.e., which
+   *   quantities can be returned in calls to geod_position() and
+   *   geod_genposition().
+   *
+   * \e g must have been initialized with a call to geod_init().  \e lat1
+   * should be in the range [&minus;90&deg;, 90&deg;]; \e azi1 should be in the
+   * range [&minus;540&deg;, 540&deg;).
+   *
+   * The geod_mask values are
+   * - \e caps |= GEOD_LATITUDE for the latitude \e lat2; this is
+   *   added automatically
+   * - \e caps |= GEOD_LONGITUDE for the latitude \e lon2
+   * - \e caps |= GEOD_AZIMUTH for the latitude \e azi2; this is
+   *   added automatically
+   * - \e caps |= GEOD_DISTANCE for the distance \e s12
+   * - \e caps |= GEOD_REDUCEDLENGTH for the reduced length \e m12
+   * - \e caps |= GEOD_GEODESICSCALE for the geodesic scales \e M12
+   *   and \e M21
+   * - \e caps |= GEOD_AREA for the area \e S12
+   * - \e caps |= GEOD_DISTANCE_IN permits the length of the
+   *   geodesic to be given in terms of \e s12; without this capability the
+   *   length can only be specified in terms of arc length.
+   * .
+   * A value of \e caps = 0 is treated as GEOD_LATITUDE | GEOD_LONGITUDE |
+   * GEOD_AZIMUTH | GEOD_DISTANCE_IN (to support the solution of the "standard"
+   * direct problem).
+   **********************************************************************/
+  void geod_lineinit(struct geod_geodesicline* l,
+                     const struct geod_geodesic* g,
+                     double lat1, double lon1, double azi1, unsigned caps);
+
+  /**
+   * Solve the direct geodesic problem.
+   *
+   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] lat1 latitude of point 1 (degrees).
+   * @param[in] lon1 longitude of point 1 (degrees).
+   * @param[in] azi1 azimuth at point 1 (degrees).
+   * @param[in] s12 distance between point 1 and point 2 (meters); it can be
+   *   negative.
+   * @param[out] plat2 pointer to the latitude of point 2 (degrees).
+   * @param[out] plon2 pointer to the longitude of point 2 (degrees).
+   * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
+   *
+   * \e g must have been initialized with a call to geod_init().  \e lat1
+   * should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e azi1
+   * should be in the range [&minus;540&deg;, 540&deg;).  The values of \e lon2
+   * and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;).  Any of
+   * the "return" arguments plat2, etc., may be replaced by 0, if you do not
+   * need some quantities computed.
+   *
+   * If either point is at a pole, the azimuth is defined by keeping the
+   * longitude fixed and writing \e lat = &plusmn;(90&deg; &minus; &epsilon;)
+   * and taking the limit &epsilon; &rarr; 0+.  An arc length greater that
+   * 180&deg; signifies a geodesic which is not a shortest path.  (For a
+   * prolate ellipsoid, an additional condition is necessary for a shortest
+   * path: the longitudinal extent must not exceed of 180&deg;.)
+   *
+   * Example, determine the point 10000 km NE of JFK:
+   @code
+   struct geod_geodesic g;
+   double lat, lon;
+   geod_init(&g, 6378137, 1/298.257223563);
+   geod_direct(&g, 40.64, -73.78, 45.0, 10e6, &lat, &lon, 0);
+   printf("%.5f %.5f\n", lat, lon);
+   @endcode
+   **********************************************************************/
+  void geod_direct(const struct geod_geodesic* g,
+                   double lat1, double lon1, double azi1, double s12,
+                   double* plat2, double* plon2, double* pazi2);
+
+  /**
+   * Solve the inverse geodesic problem.
+   *
+   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] lat1 latitude of point 1 (degrees).
+   * @param[in] lon1 longitude of point 1 (degrees).
+   * @param[in] lat2 latitude of point 2 (degrees).
+   * @param[in] lon2 longitude of point 2 (degrees).
+   * @param[out] ps12 pointer to the distance between point 1 and point 2
+   *   (meters).
+   * @param[out] pazi1 pointer to the azimuth at point 1 (degrees).
+   * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
+   *
+   * \e g must have been initialized with a call to geod_init().  \e lat1
+   * and \e lat2 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and
+   * \e lon2 should be in the range [&minus;540&deg;, 540&deg;).  The values of
+   * \e azi1 and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;).
+   * Any of the "return" arguments ps12, etc., may be replaced by 0, if you do
+   * not need some quantities computed.
+   *
+   * If either point is at a pole, the azimuth is defined by keeping the
+   * longitude fixed and writing \e lat = &plusmn;(90&deg; &minus; &epsilon;)
+   * and taking the limit &epsilon; &rarr; 0+.
+   *
+   * The solution to the inverse problem is found using Newton's method.  If
+   * this fails to converge (this is very unlikely in geodetic applications
+   * but does occur for very eccentric ellipsoids), then the bisection method
+   * is used to refine the solution.
+   *
+   * Example, determine the distance between JFK and Singapore Changi Airport:
+   @code
+   struct geod_geodesic g;
+   double s12;
+   geod_init(&g, 6378137, 1/298.257223563);
+   geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, 0, 0);
+   printf("%.3f\n", s12);
+   @endcode
+   **********************************************************************/
+  void geod_inverse(const struct geod_geodesic* g,
                     double lat1, double lon1, double lat2, double lon2,
-                    double* ps12, double* pazi1, double* pazi2,
-                    double* pm12, double* pM12, double* pM21, double* pS12);
-  double GenPosition(const struct GeodesicLine* l,
-                     int arcmode, double s12_a12,
-                     double* plat2, double* plon2, double* pazi2,
-                     double* ps12, double* pm12,
-                     double* pM12, double* pM21,
-                     double* pS12);
-
-  enum mask {
-    NONE          = 0U,
-    LATITUDE      = 1U<<7  | 0U,
-    LONGITUDE     = 1U<<8  | 1U<<3,
-    AZIMUTH       = 1U<<9  | 0U,
-    DISTANCE      = 1U<<10 | 1U<<0,
-    DISTANCE_IN   = 1U<<11 | 1U<<0 | 1U<<1,
-    REDUCEDLENGTH = 1U<<12 | 1U<<0 | 1U<<2,
-    GEODESICSCALE = 1U<<13 | 1U<<0 | 1U<<2,
-    AREA          = 1U<<14 | 1U<<4,
-    ALL           = 0x7F80U| 0x1FU
+                    double* ps12, double* pazi1, double* pazi2);
+
+  /**
+   * Compute the position along a geod_geodesicline.
+   *
+   * @param[in] l the geod_geodesicline object specifying the geodesic line.
+   * @param[in] s12 distance between point 1 and point 2 (meters); it can be
+   *   negative.
+   * @param[out] plat2 pointer to the latitude of point 2 (degrees).
+   * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires
+   *   that the \e was initialized with \e caps |= GEOD_LONGITUDE.
+   * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
+   *
+   * \e l must have been initialized with a call to geod_lineinit() with \e
+   * caps |= GEOD_DISTANCE_IN.  The values of \e lon2 and \e azi2 returned are
+   * in the range [&minus;180&deg;, 180&deg;).  Any of the "return" arguments
+   * plat2, etc., may be replaced by 0, if you do not need some quantities
+   * computed.
+   *
+   * Example, compute way points between JFK and Singapore Changi Airport
+   * the "obvious" way using geod_direct():
+   @code
+   struct geod_geodesic g;
+   double s12, azi1, lat[101],lon[101];
+   int i;
+   geod_init(&g, 6378137, 1/298.257223563);
+   geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0);
+   for (i = 0; i < 101; ++i) {
+     geod_direct(&g, 40.64, -73.78, azi1, i * s12 * 0.01, lat + i, lon + i, 0);
+     printf("%.5f %.5f\n", lat[i], lon[i]);
+   }
+   @endcode
+   * A faster way using geod_position():
+   @code
+   struct geod_geodesic g;
+   struct geod_geodesicline l;
+   double s12, azi1, lat[101],lon[101];
+   int i;
+   geod_init(&g, 6378137, 1/298.257223563);
+   geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0);
+   geod_lineinit(&l, &g, 40.64, -73.78, azi1, 0);
+   for (i = 0; i < 101; ++i) {
+     geod_position(&l, i * s12 * 0.01, lat + i, lon + i, 0);
+     printf("%.5f %.5f\n", lat[i], lon[i]);
+   }
+   @endcode
+   **********************************************************************/
+  void geod_position(const struct geod_geodesicline* l, double s12,
+                     double* plat2, double* plon2, double* pazi2);
+
+  /**
+   * The general direct geodesic problem.
+   *
+   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] lat1 latitude of point 1 (degrees).
+   * @param[in] lon1 longitude of point 1 (degrees).
+   * @param[in] azi1 azimuth at point 1 (degrees).
+   * @param[in] arcmode flag determining the meaning of the \e
+   *   s12_a12.
+   * @param[in] s12_a12 if \e arcmode is 0, this is the distance between
+   *   point 1 and point 2 (meters); otherwise it is the arc length between
+   *   point 1 and point 2 (degrees); it can be negative.
+   * @param[out] plat2 pointer to the latitude of point 2 (degrees).
+   * @param[out] plon2 pointer to the longitude of point 2 (degrees).
+   * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
+   * @param[out] ps12 pointer to the distance between point 1 and point 2
+   *   (meters).
+   * @param[out] pm12 pointer to the reduced length of geodesic (meters).
+   * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
+   *   point 1 (dimensionless).
+   * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
+   *   point 2 (dimensionless).
+   * @param[out] pS12 pointer to the area under the geodesic
+   *   (meters<sup>2</sup>).
+   * @return \e a12 arc length of between point 1 and point 2 (degrees).
+   *
+   * \e g must have been initialized with a call to geod_init().  \e lat1
+   * should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e azi1
+   * should be in the range [&minus;540&deg;, 540&deg;).  The function value \e
+   * a12 equals \e s12_a12 is \e arcmode is non-zero.  Any of the "return"
+   * arguments plat2, etc., may be replaced by 0, if you do not need some
+   * quantities computed.
+   **********************************************************************/
+  double geod_gendirect(const struct geod_geodesic* g,
+                        double lat1, double lon1, double azi1,
+                        int arcmode, double s12_a12,
+                        double* plat2, double* plon2, double* pazi2,
+                        double* ps12, double* pm12, double* pM12, double* pM21,
+                        double* pS12);
+
+  /**
+   * The general inverse geodesic calculation.
+   *
+   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] lat1 latitude of point 1 (degrees).
+   * @param[in] lon1 longitude of point 1 (degrees).
+   * @param[in] lat2 latitude of point 2 (degrees).
+   * @param[in] lon2 longitude of point 2 (degrees).
+   * @param[out] ps12 pointer to the distance between point 1 and point 2
+   *  (meters).
+   * @param[out] pazi1 pointer to the azimuth at point 1 (degrees).
+   * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
+   * @param[out] pm12 pointer to the reduced length of geodesic (meters).
+   * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
+   *   point 1 (dimensionless).
+   * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
+   *   point 2 (dimensionless).
+   * @param[out] pS12 pointer to the area under the geodesic
+   *   (meters<sup>2</sup>).
+   * @return \e a12 arc length of between point 1 and point 2 (degrees).
+   *
+   * \e g must have been initialized with a call to geod_init().  \e lat1
+   * and \e lat2 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and
+   * \e lon2 should be in the range [&minus;540&deg;, 540&deg;).  Any of the
+   * "return" arguments ps12, etc., may be replaced by 0, if you do not need
+   * some quantities computed.
+   **********************************************************************/
+  double geod_geninverse(const struct geod_geodesic* g,
+                         double lat1, double lon1, double lat2, double lon2,
+                         double* ps12, double* pazi1, double* pazi2,
+                         double* pm12, double* pM12, double* pM21,
+                         double* pS12);
+
+  /**
+   * The general position function.
+   *
+   * @param[in] l the geod_geodesicline object specifying the geodesic line.
+   * @param[in] arcmode flag determining the meaning of the second parameter;
+   *   if arcmode is 0, then \e l must have been initialized with \e caps |=
+   *   GEOD_DISTANCE_IN.
+   * @param[in] s12_a12 if \e arcmode is 0, this is the distance between
+   *   point 1 and point 2 (meters); otherwise it is the arc length between
+   *   point 1 and point 2 (degrees); it can be negative.
+   * @param[out] plat2 pointer to the latitude of point 2 (degrees).
+   * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires
+   *   that \e l was initialized with \e caps |= GEOD_LONGITUDE.
+   * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
+   * @param[out] ps12 pointer to the distance between point 1 and point 2
+   *   (meters); requires that \e l was initialized with \e caps |=
+   *   GEOD_DISTANCE.
+   * @param[out] pm12 pointer to the reduced length of geodesic (meters);
+   *   requires that \e l was initialized with \e caps |= GEOD_REDUCEDLENGTH.
+   * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
+   *   point 1 (dimensionless); requires that \e l was initialized with \e caps
+   *   |= GEOD_GEODESICSCALE.
+   * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
+   *   point 2 (dimensionless); requires that \e l was initialized with \e caps
+   *   |= GEOD_GEODESICSCALE.
+   * @param[out] pS12 pointer to the area under the geodesic
+   *   (meters<sup>2</sup>); requires that \e l was initialized with \e caps |=
+   *   GEOD_AREA.
+
+   * @return \e a12 arc length of between point 1 and point 2 (degrees).
+   *
+   * \e l must have been initialized with a call to geod_lineinit() with \e
+   * caps |= GEOD_DISTANCE_IN.  The values of \e lon2 and \e azi2 returned are
+   * in the range [&minus;180&deg;, 180&deg;).  Any of the "return" arguments
+   * plat2, etc., may be replaced by 0, if you do not need some quantities
+   * computed.  Requesting a value which \e l is not capable of computing is
+   * not an error; the corresponding argument will not be altered.
+   *
+   * Example, computer way points between JFK and Singapore Changi Airport
+   * using geod_genposition().  In this example, the points are evenly space in
+   * arc length (and so only approximately equally space in distance).  This is
+   * faster than using geod_position() would be appropriate if drawing the path
+   * on a map.
+   @code
+   struct geod_geodesic g;
+   struct geod_geodesicline l;
+   double a12, azi1, lat[101],lon[101];
+   int i;
+   geod_init(&g, 6378137, 1/298.257223563);
+   a12 = geod_geninverse(&g, 40.64, -73.78, 1.36, 103.99,
+                         0, &azi1, 0, 0, 0, 0, 0);
+   geod_lineinit(&l, &g, 40.64, -73.78, azi1, GEOD_LATITUDE | GEOD_LONGITUDE);
+   for (i = 0; i < 101; ++i) {
+     geod_genposition(&l, 1, i * a12 * 0.01,
+                      lat + i, lon + i, 0, 0, 0, 0, 0, 0);
+     printf("%.5f %.5f\n", lat[i], lon[i]);
+   }
+   @endcode
+   **********************************************************************/
+  double geod_genposition(const struct geod_geodesicline* l,
+                          int arcmode, double s12_a12,
+                          double* plat2, double* plon2, double* pazi2,
+                          double* ps12, double* pm12,
+                          double* pM12, double* pM21,
+                          double* pS12);
+
+  /**
+   * The area of a geodesic polygon.
+   *
+   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] lats an array of latitudes of the polygon vertices (degrees).
+   * @param[in] lons an array of longitudes of the polygon vertices (degrees).
+   * @param[in] n the number of vertices.
+   * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>).
+   * @param[out] pP pointer to the perimeter of the polygon (meters).
+   *
+   * \e lats should be in the range [&minus;90&deg;, 90&deg;]; \e lons should
+   * be in the range [&minus;540&deg;, 540&deg;).
+   *
+   * Only simple polygons (which are not self-intersecting) are allowed.
+   * There's no need to "close" the polygon by repeating the first vertex.  The
+   * area returned is signed with counter-clockwise traversal being treated as
+   * positive.
+   *
+   * Example, compute the area of Antarctic:
+   @code
+   double
+     lats[] = {-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
+               -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7},
+     lons[] = {-74, -102, -102, -131, -163, 163, 172, 140, 113,
+                88, 59, 25, -4, -14, -33, -46, -61};
+   struct geod_geodesic g;
+   double A, P;
+   geod_init(&g, 6378137, 1/298.257223563);
+   geod_polygonarea(&g, lats, lons, (sizeof lats) / (sizeof lats[0]), &A, &P);
+   printf("%.0f %.2f\n", A, P);
+   @endcode
+   **********************************************************************/
+  void geod_polygonarea(const struct geod_geodesic* g,
+                        double lats[], double lons[], int n,
+                        double* pA, double* pP);
+
+  /**
+   * mask values for geod_geodesicline capabilities.
+   **********************************************************************/
+  enum geod_mask {
+    GEOD_NONE         = 0U,
+    GEOD_LATITUDE     = 1U<<7  | 0U,            /**< Calculate latitude */
+    GEOD_LONGITUDE    = 1U<<8  | 1U<<3,         /**< Calculate longitude */
+    GEOD_AZIMUTH      = 1U<<9  | 0U,            /**< Calculate azimuth */
+    GEOD_DISTANCE     = 1U<<10 | 1U<<0,         /**< Calculate distance */
+    GEOD_DISTANCE_IN  = 1U<<11 | 1U<<0 | 1U<<1, /**< Allow distance as input  */
+    GEOD_REDUCEDLENGTH= 1U<<12 | 1U<<0 | 1U<<2, /**< Calculate reduced length */
+    GEOD_GEODESICSCALE= 1U<<13 | 1U<<0 | 1U<<2, /**< Calculate geodesic scale */
+    GEOD_AREA         = 1U<<14 | 1U<<4,         /**< Calculate reduced length */
+    GEOD_ALL          = 0x7F80U| 0x1FU          /**< Calculate everything */
   };
 
 #if defined(__cplusplus)