Fix NaN handling by geod_inverse and incorrect area computation by geod_polygon_addedge if the longitude extent of the added edge exceeds 180 degrees (patch by Charles Karney, #251, #253)

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

2 files changed, 120 insertions, 73 deletions

diff --git a/src/geodesic.c b/src/geodesic.c
index bd9fc960..fd0214c7 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -13,12 +13,12 @@
  *   C. F. F. Karney,
  *   Algorithms for geodesics,
  *   J. Geodesy <b>87</b>, 43--55 (2013);
- *   http://dx.doi.org/10.1007/s00190-012-0578-z
+ *   https://dx.doi.org/10.1007/s00190-012-0578-z
  *   Addenda: http://geographiclib.sf.net/geod-addenda.html
  *
  * See the comments in geodesic.h for documentation.
  *
- * Copyright (c) Charles Karney (2012-2013) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2014) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * http://geographiclib.sourceforge.net/
  */
@@ -222,6 +222,7 @@ static void C1pf(real eps, real c[]);
 static real A2m1f(real eps);
 static void C2f(real eps, real c[]);
 static int transit(real lon1, real lon2);
+static int transitdirect(real lon1, real lon2);
 static void accini(real s[]);
 static void acccopy(const real s[], real t[]);
 static void accadd(real s[], real y);
@@ -271,18 +272,16 @@ void geod_lineinit(struct geod_geodesicline* l,
   l->caps = (caps ? caps : GEOD_DISTANCE_IN | GEOD_LONGITUDE) |
     GEOD_LATITUDE | GEOD_AZIMUTH; /* Always allow latitude and azimuth */
 
-  /* Guard against underflow in salp0 */
-  azi1 = AngRound(AngNormalize(azi1));
-  lon1 = AngNormalize(lon1);
   l->lat1 = lat1;
   l->lon1 = lon1;
-  l->azi1 = azi1;
+  /* Guard against underflow in salp0 */
+  l->azi1 = AngRound(AngNormalize(azi1));
   /* alp1 is in [0, pi] */
-  alp1 = azi1 * degree;
+  alp1 = l->azi1 * degree;
   /* Enforce sin(pi) == 0 and cos(pi/2) == 0.  Better to face the ensuing
    * problems directly than to skirt them. */
-  l->salp1 =      azi1  == -180 ? 0 : sin(alp1);
-  l->calp1 = fabs(azi1) ==   90 ? 0 : cos(alp1);
+  l->salp1 =      l->azi1  == -180 ? 0 : sin(alp1);
+  l->calp1 = fabs(l->azi1) ==   90 ? 0 : cos(alp1);
   phi = lat1 * degree;
   /* Ensure cbet1 = +epsilon at poles */
   sbet1 = l->f1 * sin(phi);
@@ -349,7 +348,7 @@ void geod_lineinit(struct geod_geodesicline* l,
 }
 
 real geod_genposition(const struct geod_geodesicline* l,
-                      boolx arcmode, real s12_a12,
+                      unsigned flags, real s12_a12,
                       real* plat2, real* plon2, real* pazi2,
                       real* ps12, real* pm12,
                       real* pM12, real* pM21,
@@ -371,11 +370,11 @@ real geod_genposition(const struct geod_geodesicline* l,
 
   outmask &= l->caps & OUT_ALL;
   if (!( TRUE /*Init()*/ &&
-         (arcmode || (l->caps & GEOD_DISTANCE_IN & OUT_ALL)) ))
+         (flags & GEOD_ARCMODE || (l->caps & GEOD_DISTANCE_IN & OUT_ALL)) ))
     /* Uninitialized or impossible distance calculation requested */
     return NaN;
 
-  if (arcmode) {
+  if (flags & GEOD_ARCMODE) {
     real s12a;
     /* Interpret s12_a12 as spherical arc length */
     sig12 = s12_a12 * degree;
@@ -435,7 +434,7 @@ real geod_genposition(const struct geod_geodesicline* l,
   csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;
   dn2 = sqrt(1 + l->k2 * sq(ssig2));
   if (outmask & (GEOD_DISTANCE | GEOD_REDUCEDLENGTH | GEOD_GEODESICSCALE)) {
-    if (arcmode || fabs(l->f) > 0.01)
+    if (flags & GEOD_ARCMODE || fabs(l->f) > 0.01)
       B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
     AB1 = (1 + l->A1m1) * (B12 - l->B11);
   }
@@ -446,26 +445,29 @@ real geod_genposition(const struct geod_geodesicline* l,
   if (cbet2 == 0)
     /* I.e., salp0 = 0, csig2 = 0.  Break the degeneracy in this case */
     cbet2 = csig2 = tiny;
-  /* tan(omg2) = sin(alp0) * tan(sig2) */
-  somg2 = l->salp0 * ssig2; comg2 = csig2;  /* No need to normalize */
   /* tan(alp0) = cos(sig2)*tan(alp2) */
   salp2 = l->salp0; calp2 = l->calp0 * csig2; /* No need to normalize */
-  /* omg12 = omg2 - omg1 */
-  omg12 = atan2(somg2 * l->comg1 - comg2 * l->somg1,
-                comg2 * l->comg1 + somg2 * l->somg1);
 
   if (outmask & GEOD_DISTANCE)
-    s12 = arcmode ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12;
+    s12 = flags & GEOD_ARCMODE ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12;
 
   if (outmask & GEOD_LONGITUDE) {
+    /* tan(omg2) = sin(alp0) * tan(sig2) */
+    somg2 = l->salp0 * ssig2; comg2 = csig2;  /* No need to normalize */
+    /* omg12 = omg2 - omg1 */
+    omg12 = flags & GEOD_LONG_NOWRAP ? sig12
+      - (atan2(ssig2, csig2) - atan2(l->ssig1, l->csig1))
+      + (atan2(somg2, comg2) - atan2(l->somg1, l->comg1))
+      : atan2(somg2 * l->comg1 - comg2 * l->somg1,
+              comg2 * l->comg1 + somg2 * l->somg1);
     lam12 = omg12 + l->A3c *
       ( sig12 + (SinCosSeries(TRUE, ssig2, csig2, l->C3a, nC3-1)
                  - l->B31));
     lon12 = lam12 / degree;
     /* Use AngNormalize2 because longitude might have wrapped multiple
      * times. */
-    lon12 = AngNormalize2(lon12);
-    lon2 = AngNormalize(l->lon1 + lon12);
+    lon2 = flags & GEOD_LONG_NOWRAP ? l->lon1 + lon12 :
+      AngNormalize(AngNormalize(l->lon1) + AngNormalize2(lon12));
   }
 
   if (outmask & GEOD_LATITUDE)
@@ -542,7 +544,7 @@ real geod_genposition(const struct geod_geodesicline* l,
   if (outmask & GEOD_AREA)
     *pS12 = S12;
 
-  return arcmode ? s12_a12 : sig12 / degree;
+  return flags & GEOD_ARCMODE ? s12_a12 : sig12 / degree;
 }
 
 void geod_position(const struct geod_geodesicline* l, real s12,
@@ -552,7 +554,7 @@ void geod_position(const struct geod_geodesicline* l, real s12,
 
 real geod_gendirect(const struct geod_geodesic* g,
                     real lat1, real lon1, real azi1,
-                    boolx arcmode, real s12_a12,
+                    unsigned flags, real s12_a12,
                     real* plat2, real* plon2, real* pazi2,
                     real* ps12, real* pm12, real* pM12, real* pM21,
                     real* pS12) {
@@ -568,8 +570,9 @@ real geod_gendirect(const struct geod_geodesic* g,
 
   geod_lineinit(&l, g, lat1, lon1, azi1,
                 /* Automatically supply GEOD_DISTANCE_IN if necessary */
-                outmask | (arcmode ? GEOD_NONE : GEOD_DISTANCE_IN));
-  return geod_genposition(&l, arcmode, s12_a12,
+                outmask |
+                (flags & GEOD_ARCMODE ? GEOD_NONE : GEOD_DISTANCE_IN));
+  return geod_genposition(&l, flags, s12_a12,
                           plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12);
 }
 
@@ -577,7 +580,7 @@ void geod_direct(const struct geod_geodesic* g,
                  real lat1, real lon1, real azi1,
                  real s12,
                  real* plat2, real* plon2, real* pazi2) {
-  geod_gendirect(g, lat1, lon1, azi1, FALSE, s12, plat2, plon2, pazi2,
+  geod_gendirect(g, lat1, lon1, azi1, GEOD_NOFLAGS, s12, plat2, plon2, pazi2,
                  0, 0, 0, 0, 0);
 }
 
@@ -700,7 +703,7 @@ real geod_geninverse(const struct geod_geodesic* g,
     /* Add the check for sig12 since zero length geodesics might yield m12 <
      * 0.  Test case was
      *
-     *    echo 20.001 0 20.001 0 | Geod -i
+     *    echo 20.001 0 20.001 0 | GeodSolve -i
      *
      * In fact, we will have sig12 > pi/2 for meridional geodesic which is
      * not a shortest path. */
@@ -1214,7 +1217,8 @@ real InverseStart(const struct geod_geodesic* g,
       calp1 = sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
     }
   }
-  if (salp1 > 0)              /* Sanity check on starting guess */
+  /* Sanity check on starting guess.  Backwards check allows NaN through. */
+  if (!(salp1 <= 0))
     SinCosNorm(&salp1, &calp1);
   else {
     salp1 = 1; calp1 = 0;
@@ -1516,6 +1520,13 @@ int transit(real lon1, real lon2) {
     (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0);
 }
 
+int transitdirect(real lon1, real lon2) {
+  lon1 = fmod(lon1, (real)(720));
+  lon2 = fmod(lon2, (real)(720));
+  return ( ((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) -
+           ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1) );
+}
+
 void accini(real s[]) {
   /* Initialize an accumulator; this is an array with two elements. */
   s[0] = s[1] = 0;
@@ -1583,13 +1594,13 @@ void geod_polygon_addedge(const struct geod_geodesic* g,
                           real azi, real s) {
   if (p->num) {                 /* Do nothing is num is zero */
     real lat, lon, S12;
-    geod_gendirect(g, p->lat, p->lon, azi, FALSE, s,
+    geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_NOWRAP, s,
                    &lat, &lon, 0,
                    0, 0, 0, 0, p->polyline ? 0 : &S12);
     accadd(p->P, s);
     if (!p->polyline) {
       accadd(p->A, S12);
-      p->crossings += transit(p->lon, lon);
+      p->crossings += transitdirect(p->lon, lon);
     }
     p->lat = lat; p->lon = lon;
     ++p->num;
@@ -1720,11 +1731,11 @@ unsigned geod_polygon_testedge(const struct geod_geodesic* g,
   crossings = p->crossings;
   {
     real lat, lon, s12, S12;
-    geod_gendirect(g, p->lat, p->lon, azi, FALSE, s,
+    geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_NOWRAP, s,
                    &lat, &lon, 0,
                    0, 0, 0, 0, &S12);
     tempsum += S12;
-    crossings += transit(p->lon, lon);
+    crossings += transitdirect(p->lon, lon);
     geod_geninverse(g, lat,  lon, p->lat0,  p->lon0,
                     &s12, 0, 0, 0, 0, 0, &S12);
     perimeter += s12;
diff --git a/src/geodesic.h b/src/geodesic.h
index 69985009..6b2afc5a 100644
--- a/src/geodesic.h
+++ b/src/geodesic.h
@@ -4,10 +4,10 @@
  *
  * 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">
+ *   <a href="https://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">
+ *   DOI: <a href="https://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>.
@@ -75,18 +75,18 @@
  * (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
+ * - \e lat1 = &minus;\e lat2 (with neither point 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.
+ * - \e lon2 = \e lon1 &plusmn; 180&deg; (with neither point 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.)
+ *   setting [\e azi1, \e azi2] = [&minus;\e azi1, &minus;\e azi2], \e S12
+ *   = &minus;\e S12.  (This occurs when \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
@@ -103,17 +103,17 @@
  * to the member functions.  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
+ * GeographicLib::GeodesicLine, and GeographicLib::PolygonAreaT, 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.
  *
- * Copyright (c) Charles Karney (2012-2013) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2014) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * http://geographiclib.sourceforge.net/
  *
  * This library was distributed with
- * <a href="../index.html">GeographicLib</a> 1.32.
+ * <a href="../index.html">GeographicLib</a> 1.40.
  **********************************************************************/
 
 #if !defined(GEODESIC_H)
@@ -128,7 +128,7 @@
  * The minor version of the geodesic library.  (This tracks the version of
  * GeographicLib.)
  **********************************************************************/
-#define GEODESIC_VERSION_MINOR 32
+#define GEODESIC_VERSION_MINOR 40
 /**
  * The patch level of the geodesic library.  (This tracks the version of
  * GeographicLib.)
@@ -381,11 +381,13 @@ extern "C" {
    * @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[in] flags bitor'ed combination of geod_flags(); \e flags &
+   *   GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags &
+   *   GEOD_LONG_NOWRAP prevents the value of \e lon2 being wrapped into
+   *   the range [&minus;180&deg;, 180&deg;).
+   * @param[in] s12_a12 if \e flags & GEOD_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).
@@ -402,14 +404,21 @@ extern "C" {
    *
    * \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 \e plat2, etc., may be replaced by 0, if you do not need some
-   * quantities computed.
+   * should be in the range [&minus;540&deg;, 540&deg;).  The function
+   * value \e a12 equals \e s12_a12 if \e flags & GEOD_ARCMODE.  Any of the
+   * "return" arguments \e plat2, etc., may be replaced by 0, if you do not
+   * need some quantities computed.
+   *
+   * With \e flags & GEOD_LONG_NOWRAP bit set, the quantity \e lon2 &minus;
+   * \e lon1 indicates how many times the geodesic wrapped around the
+   * ellipsoid.  Because \e lon2 might be outside the normal allowed range
+   * for longitudes, [&minus;540&deg;, 540&deg;), be sure to normalize it,
+   * e.g., with fmod(\e lon2, 360.0) before using it in subsequent
+   * calculations
    **********************************************************************/
   double geod_gendirect(const struct geod_geodesic* g,
                         double lat1, double lon1, double azi1,
-                        int arcmode, double s12_a12,
+                        unsigned flags, double s12_a12,
                         double* plat2, double* plon2, double* pazi2,
                         double* ps12, double* pm12, double* pM12, double* pM21,
                         double* pS12);
@@ -453,12 +462,16 @@ extern "C" {
    *
    * @param[in] l a pointer to 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 |=
+   * @param[in] flags bitor'ed combination of geod_flags(); \e flags &
+   *   GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags &
+   *   GEOD_LONG_NOWRAP prevents the value of \e lon2 being wrapped into
+   *   the range [&minus;180&deg;, 180&deg;); if \e flags & GEOD_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[in] s12_a12 if \e flags & GEOD_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.
@@ -480,11 +493,18 @@ extern "C" {
    * @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
-   * \e 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.
+   * caps |= GEOD_DISTANCE_IN.  The value \e azi2 returned is in the range
+   * [&minus;180&deg;, 180&deg;).  Any of the "return" arguments \e 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.
+   *
+   * With \e flags & GEOD_LONG_NOWRAP bit set, the quantity \e lon2 &minus;
+   * \e lon1 indicates how many times the geodesic wrapped around the
+   * ellipsoid.  Because \e lon2 might be outside the normal allowed range
+   * for longitudes, [&minus;540&deg;, 540&deg;), be sure to normalize it,
+   * e.g., with fmod(\e lon2, 360.0) before using it in subsequent
+   * calculations
    *
    * Example, compute way points between JFK and Singapore Changi Airport
    * using geod_genposition().  In this example, the points are evenly space in
@@ -494,7 +514,7 @@ extern "C" {
    @code
    struct geod_geodesic g;
    struct geod_geodesicline l;
-   double a12, azi1, lat[101],lon[101];
+   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,
@@ -508,7 +528,7 @@ extern "C" {
    @endcode
    **********************************************************************/
   double geod_genposition(const struct geod_geodesicline* l,
-                          int arcmode, double s12_a12,
+                          unsigned flags, double s12_a12,
                           double* plat2, double* plon2, double* pazi2,
                           double* ps12, double* pm12,
                           double* pM12, double* pM21,
@@ -526,6 +546,10 @@ extern "C" {
    * polylinep is non-zero, then the vertices and edges define a polyline and
    * only the perimeter is returned by geod_polygon_compute().
    *
+   * The area and perimeter are accumulated at two times the standard floating
+   * point precision to guard against the loss of accuracy with many-sided
+   * polygons.  At any point you can ask for the perimeter and area so far.
+   *
    * An example of the use of this function is given in the documentation for
    * geod_polygon_compute().
    **********************************************************************/
@@ -592,10 +616,12 @@ extern "C" {
    *   polyline (meters).
    * @return the number of points.
    *
-   * Only simple polygons (which are not self-intersecting) are allowed.
-   * There's no need to "close" the polygon by repeating the first vertex.  Set
-   * \e pA or \e pP to zero, if you do not want the corresponding quantity
-   * returned.
+   * The area and perimeter are accumulated at two times the standard floating
+   * point precision to guard against the loss of accuracy with many-sided
+   * polygons.  Only simple polygons (which are not self-intersecting) are
+   * allowed.  There's no need to "close" the polygon by repeating the first
+   * vertex.  Set \e pA or \e pP to zero, if you do not want the corresponding
+   * quantity returned.
    *
    * Example, compute the perimeter and area of the geodesic triangle with
    * vertices (0&deg;N,0&deg;E), (0&deg;N,90&deg;E), (90&deg;N,0&deg;E).
@@ -704,7 +730,7 @@ extern "C" {
    * area returned is signed with counter-clockwise traversal being treated as
    * positive.
    *
-   * Example, compute the area of Antarctic:
+   * Example, compute the area of Antarctica:
    @code
    double
      lats[] = {-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
@@ -723,7 +749,7 @@ extern "C" {
                         double* pA, double* pP);
 
   /**
-   * mask values for the the \e caps argument to geod_lineinit().
+   * mask values for the \e caps argument to geod_lineinit().
    **********************************************************************/
   enum geod_mask {
     GEOD_NONE         = 0U,                     /**< Calculate nothing */
@@ -738,6 +764,16 @@ extern "C" {
     GEOD_ALL          = 0x7F80U| 0x1FU          /**< Calculate everything */
   };
 
+  /**
+   * flag values for the \e flags argument to geod_gendirect() and
+   * geod_genposition()
+   **********************************************************************/
+  enum geod_flags {
+    GEOD_NOFLAGS      = 0U,     /**< No flags */
+    GEOD_ARCMODE      = 1U<<0,  /**< Position given in terms of arc distance */
+    GEOD_LONG_NOWRAP  = 1U<<15  /**< Don't wrap longitude */
+  };
+
 #if defined(__cplusplus)
 }
 #endif