allow polygon vertices to be specified incrementally for geodesic area (#221).

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

5 files changed, 546 insertions, 105 deletions

diff --git a/ChangeLog b/ChangeLog
index 231490e3..73859021 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2013-07-12  Frank Warmerdam  <warmerdam@pobox.com>
+
+	* src/geodesic.{c,h}: allow polygon vertices to be specified 
+	incrementally for geodesic area (#221).
+
 2013-07-08  Frank Warmerdam  <warmerdam@pobox.com>
 
 	* src/PJ_calcofi.c: Add Cal Coop Ocean Fish Invest Lines/Stations
diff --git a/man/man1/geod.1 b/man/man1/geod.1
index 0a2522ea..8ff59db2 100644
--- a/man/man1/geod.1
+++ b/man/man1/geod.1
@@ -1,8 +1,8 @@
 .\" @(#)geod.1
-.nr LL 6.5i
+.nr LL 7.0i
 .ad b
 .hy 1
-.TH GEOD 1 "2013/02/27 Rel. 4.8"
+.TH GEOD 1 "2013/07/11 Rel. 4.9.0"
 .SH NAME
 geod \- direct geodesic computations
 .br
@@ -38,9 +38,10 @@ perform geodesic (\(``Great Circle\('') computations for determining
 latitude, longitude and back azimuth of a terminus point
 given a initial point latitude, longitude, azimuth and distance (direct) or
 the forward and back azimuths and distance between an initial and
-terminus point latitudes and longitudes (inverse).
+terminus point latitudes and longitudes (inverse).  The results are
+accurate to round off for |\fIf\fR| < 1/50, where \fIf\fR is flattening.
 .PP
-The following runline control parameters can appear in any order:
+The following command-line options can appear in any order:
 .TP
 .B \-I
 Specifies that the inverse geodesic computation is to be performed.
@@ -93,13 +94,13 @@ DMS numbers between 0 and 360 degrees.  Also note -f.
 .PP
 The
 .B +args
-run-line arguments are associated with geodetic parameters
+command-line options are associated with geodetic parameters
 for specifying the ellipsoidal or sphere to use.
 See
 .B proj
 documentation for full list of these parameters and controls.
 The options are processed in left to right order
-from the run line.
+from the command line.
 Reentry of an option is ignored with the first occurrence assumed to
 be the desired value.
 .PP
@@ -200,18 +201,23 @@ Note: lack of precision in the distance value compromises
 the precision of the Portland location.
 .SH SEE ALSO
 .B geodesic(3)
-.br
-The \fBGeodSolve\fR utility in \fBGeographicLib\fR,
-http://geographiclib.sf.net.  With the \fB-E\fR option, this solves the
-geodesic problems in terms of elliptic integrals; the results are
-accurate for arbitrary \fIf\fR.
-.br
+.PP
+\fBGeographicLib\fR, http://geographiclib.sf.net
+.PP
+The \fBGeodSolve\fR utility in GeographicLib.  With the \fB-E\fR option,
+this solves the geodesic problems in terms of elliptic integrals; the
+results are accurate for arbitrary \fIf\fR.
+.PP
 C. F. F. Karney, \fIAlgorithms for Geodesics\fR,
 .br
 J. Geodesy \fB87\fR, 43-55 (2013);
 .br
-DOI: http://dx.doi.org/10.1007/s00190-012-0578-z,
+DOI: http://dx.doi.org/10.1007/s00190-012-0578-z
+.br
+http://geographiclib.sf.net/geod-addenda.html
+.PP
+The \fIonline geodesic bibliography\fR,
 .br
-http://geographiclib.sf.net/geod-addenda.html.
+http://geographiclib.sf.net/geodesic-papers/biblio.html
 .SH HOME PAGE
 http://proj.osgeo.org
diff --git a/man/man3/geodesic.3 b/man/man3/geodesic.3
index e409079d..717be2b9 100644
--- a/man/man3/geodesic.3
+++ b/man/man3/geodesic.3
@@ -1,6 +1,6 @@
 .\" @(#)geodesic.3
-.nr LL 6.5i
-.TH GEODESIC 3 "2013/02/27 Rel. 4.8" 
+.nr LL 7.0i
+.TH GEODESIC 3 "2013/07/11 Rel. 4.9.0" 
 .ad b
 .hy 1
 .SH NAME
@@ -32,13 +32,13 @@ This library is a port of the geodesic routines in the C++ library,
 GeographicLib, to C.  It solves the direct and inverse geodesic problems
 on an ellipsoid of revolution.  In addition, the reduced length of a
 geodesic and the area between a geodesic and the equator can be
-computed.  The results are accurate to round off for |\fIf\fR| < 1/50.
-Note that the geodesic routines measure angles (latitudes, longitudes,
-and azimuths) in degrees, unlike the rest of the \fBproj\fR library,
-which uses radians.  The documentation for this library is included in
-geodesic.h.  A formatted version of the documentation is available at
-.br
-http://geographiclib.sf.net/html/C/index.html
+computed.  The results are accurate to round off for |\fIf\fR| < 1/50,
+where \fIf\fR is the flattening.  Note that the geodesic routines
+measure angles (latitudes, longitudes, and azimuths) in degrees, unlike
+the rest of the \fBproj\fR library, which uses radians.  The
+documentation for this library is included in geodesic.h.  A formatted
+version of the documentation is available at
+http://geographiclib.sf.net/1.32/C
 .SH EXAMPLE
 The following program reads in lines with the coordinates for two points
 in decimal degrees (\fIlat1\fR, \fIlon1\fR, \fIlat2\fR, \fIlon2\fR) and
@@ -70,18 +70,28 @@ int main() {
 .SH LIBRARY
 libproj.a \- library of projections and support procedures
 .SH SEE ALSO
-.B geod(1)
-.br
-The \fBGeodesicExact\fR class in \fBGeographicLib\fR,
-http://geographiclib.sf.net.  This solves the geodesic problems in terms
-of elliptic integrals; the results are accurate for arbitrary \fIf\fR.
+Full online documentation for \fBgeodesic(3)\fR,
 .br
+http://geographiclib.sf.net/1.32/C
+.PP
+.B geod(1)
+.PP
+\fBGeographicLib\fR, http://geographiclib.sf.net
+.PP
+The \fBGeodesicExact\fR class in GeographicLib solves the geodesic
+problems in terms of elliptic integrals; the results are accurate for
+arbitrary \fIf\fR.
+.PP
 C. F. F. Karney, \fIAlgorithms for Geodesics\fR,
 .br
 J. Geodesy \fB87\fR, 43-55 (2013);
 .br
-DOI: http://dx.doi.org/10.1007/s00190-012-0578-z,
+DOI: http://dx.doi.org/10.1007/s00190-012-0578-z
+.br
+http://geographiclib.sf.net/geod-addenda.html
+.PP
+The \fIonline geodesic bibliography\fR,
 .br
-http://geographiclib.sf.net/geod-addenda.html.
+http://geographiclib.sf.net/geodesic-papers/biblio.html
 .SH HOME PAGE
 http://proj.osgeo.org
diff --git a/src/geodesic.c b/src/geodesic.c
index 957f2144..bd9fc960 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -1,8 +1,12 @@
 /**
  * \file geodesic.c
  * \brief Implementation of the geodesic routines in C
+ *
+ * For the full documentation see geodesic.h.
  **********************************************************************/
 
+/** @cond SKIP */
+
 /*
  * This is a C implementation of the geodesic algorithms described in
  *
@@ -20,8 +24,6 @@
  */
 
 #include "geodesic.h"
-
-/** @cond SKIP */
 #include <math.h>
 
 #define GEOGRAPHICLIB_GEODESIC_ORDER 6
@@ -221,9 +223,10 @@ static real A2m1f(real eps);
 static void C2f(real eps, real c[]);
 static int transit(real lon1, real lon2);
 static void accini(real s[]);
+static void acccopy(const real s[], real t[]);
 static void accadd(real s[], real y);
-
-/** @endcond */
+static real accsum(const real s[], real y);
+static void accneg(real s[]);
 
 void geod_init(struct geod_geodesic* g, real a, real f) {
   if (!init) Init();
@@ -925,16 +928,12 @@ real geod_geninverse(const struct geod_geodesic* g,
 }
 
 void geod_inverse(const struct geod_geodesic* g,
-                  double lat1, double lon1, double lat2, double lon2,
-                  double* ps12, double* pazi1, double* pazi2) {
+                  real lat1, real lon1, real lat2, real lon2,
+                  real* ps12, real* pazi1, real* pazi2) {
   geod_geninverse(g, lat1, lon1, lat2, lon2, ps12, pazi1, pazi2, 0, 0, 0, 0);
 }
 
-/** @cond SKIP */
-
-real SinCosSeries(boolx sinp,
-                  real sinx, real cosx,
-                  const real c[], int n) {
+real SinCosSeries(boolx sinp, real sinx, real cosx, const real c[], int n) {
   /* Evaluate
    * y = sinp ? sum(c[i] * sin( 2*i    * x), i, 1, n) :
    *            sum(c[i] * cos((2*i+1) * x), i, 0, n-1)
@@ -1522,6 +1521,11 @@ void accini(real s[]) {
   s[0] = s[1] = 0;
 }
 
+void acccopy(const real s[], real t[]) {
+  /* Copy an accumulator; t = s. */
+  t[0] = s[0]; t[1] = s[1];
+}
+
 void accadd(real s[], real y) {
   /* Add y to an accumulator. */
   real u, z = sumx(y, s[1], &u);
@@ -1532,30 +1536,235 @@ void accadd(real s[], real y) {
     s[1] = s[1] + u;
 }
 
-/** @endcond */
+real accsum(const real s[], real y) {
+  /* Return accumulator + y (but don't add to accumulator). */
+  real t[2];
+  acccopy(s, t);
+  accadd(t, y);
+  return t[0];
+}
 
-void geod_polygonarea(const struct geod_geodesic* g,
-                      real lats[], real lons[], int n,
-                      real* pA, real* pP) {
-  int i, crossings = 0;
-  real area0 = 4 * pi * g->c2, A[2], P[2];
-  accini(A); accini(P);
-  for (i = 0; i < n; ++i) {
+void accneg(real s[]) {
+  /* Negate an accumulator. */
+  s[0] = -s[0]; s[1] = -s[1];
+}
+
+void geod_polygon_init(struct geod_polygon* p, boolx polylinep) {
+  p->lat0 = p->lon0 = p->lat = p->lon = NaN;
+  p->polyline = (polylinep != 0);
+  accini(p->P);
+  accini(p->A);
+  p->num = p->crossings = 0;
+}
+
+void geod_polygon_addpoint(const struct geod_geodesic* g,
+                           struct geod_polygon* p,
+                           real lat, real lon) {
+  lon = AngNormalize(lon);
+  if (p->num == 0) {
+    p->lat0 = p->lat = lat;
+    p->lon0 = p->lon = lon;
+  } else {
+    real s12, S12;
+    geod_geninverse(g, p->lat, p->lon, lat, lon,
+                    &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12);
+    accadd(p->P, s12);
+    if (!p->polyline) {
+      accadd(p->A, S12);
+      p->crossings += transit(p->lon, lon);
+    }
+    p->lat = lat; p->lon = lon;
+  }
+  ++p->num;
+}
+
+void geod_polygon_addedge(const struct geod_geodesic* g,
+                          struct geod_polygon* p,
+                          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,
+                   &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->lat = lat; p->lon = lon;
+    ++p->num;
+  }
+}
+
+unsigned geod_polygon_compute(const struct geod_geodesic* g,
+                              const struct geod_polygon* p,
+                              boolx reverse, boolx sign,
+                              real* pA, real* pP) {
+  real s12, S12, t[2], area0;
+  int crossings;
+  if (p->num < 2) {
+    if (pP) *pP = 0;
+    if (!p->polyline && pA) *pA = 0;
+    return p->num;
+  }
+  if (p->polyline) {
+    if (pP) *pP = p->P[0];
+    return p->num;
+  }
+  geod_geninverse(g, p->lat, p->lon, p->lat0, p->lon0,
+                  &s12, 0, 0, 0, 0, 0, &S12);
+  if (pP) *pP = accsum(p->P, s12);
+  acccopy(p->A, t);
+  accadd(t, S12);
+  crossings = p->crossings + transit(p->lon, p->lon0);
+  area0 = 4 * pi * g->c2;
+  if (crossings & 1)
+    accadd(t, (t[0] < 0 ? 1 : -1) * area0/2);
+  /* area is with the clockwise sense.  If !reverse convert to
+   * counter-clockwise convention. */
+  if (!reverse)
+    accneg(t);
+  /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
+  if (sign) {
+    if (t[0] > area0/2)
+      accadd(t, -area0);
+    else if (t[0] <= -area0/2)
+      accadd(t, +area0);
+  } else {
+    if (t[0] >= area0)
+      accadd(t, -area0);
+    else if (t[0] < 0)
+      accadd(t, +area0);
+  }
+  if (pA) *pA = 0 + t[0];
+  return p->num;
+}
+
+unsigned geod_polygon_testpoint(const struct geod_geodesic* g,
+                                const struct geod_polygon* p,
+                                real lat, real lon,
+                                boolx reverse, boolx sign,
+                                real* pA, real* pP) {
+  real perimeter, tempsum, area0;
+  int crossings, i;
+  unsigned num = p->num + 1;
+  if (num == 1) {
+    if (pP) *pP = 0;
+    if (!p->polyline && pA) *pA = 0;
+    return num;
+  }
+  perimeter = p->P[0];
+  tempsum = p->polyline ? 0 : p->A[0];
+  crossings = p->crossings;
+  for (i = 0; i < (p->polyline ? 1 : 2); ++i) {
     real s12, S12;
-    geod_geninverse(g, lats[i], lons[i], lats[(i + 1) % n], lons[(i + 1) % n],
-               &s12, 0, 0, 0, 0, 0, &S12);
-    accadd(P, s12);
-    /* The minus sign is due to the counter-clockwise convention */
-    accadd(A, -S12);
-    crossings += transit(lons[i], lons[(i + 1) % n]);
+    geod_geninverse(g,
+                    i == 0 ? p->lat  : lat, i == 0 ? p->lon  : lon,
+                    i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon,
+                    &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12);
+    perimeter += s12;
+    if (!p->polyline) {
+      tempsum += S12;
+      crossings += transit(i == 0 ? p->lon  : lon,
+                           i != 0 ? p->lon0 : lon);
+    }
   }
+
+  if (pP) *pP = perimeter;
+  if (p->polyline)
+    return num;
+
+  area0 = 4 * pi * g->c2;
   if (crossings & 1)
-    accadd(A, (A[0] < 0 ? 1 : -1) * area0/2);
-  /* Put area in (-area0/2, area0/2] */
-  if (A[0] > area0/2)
-    accadd(A, -area0);
-  else if (A[0] <= -area0/2)
-    accadd(A, +area0);
-  if (pA) *pA = A[0];
-  if (pP) *pP = P[0];
+    tempsum += (tempsum < 0 ? 1 : -1) * area0/2;
+  /* area is with the clockwise sense.  If !reverse convert to
+   * counter-clockwise convention. */
+  if (!reverse)
+    tempsum *= -1;
+  /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
+  if (sign) {
+    if (tempsum > area0/2)
+      tempsum -= area0;
+    else if (tempsum <= -area0/2)
+      tempsum += area0;
+  } else {
+    if (tempsum >= area0)
+      tempsum -= area0;
+    else if (tempsum < 0)
+      tempsum += area0;
+  }
+  if (pA) *pA = 0 + tempsum;
+  return num;
 }
+
+unsigned geod_polygon_testedge(const struct geod_geodesic* g,
+                               const struct geod_polygon* p,
+                               real azi, real s,
+                               boolx reverse, boolx sign,
+                               real* pA, real* pP) {
+  real perimeter, tempsum, area0;
+  int crossings;
+  unsigned num = p->num + 1;
+  if (num == 1) {               /* we don't have a starting point! */
+    if (pP) *pP = NaN;
+    if (!p->polyline && pA) *pA = NaN;
+    return 0;
+  }
+  perimeter = p->P[0] + s;
+  if (p->polyline) {
+    if (pP) *pP = perimeter;
+    return num;
+  }
+
+  tempsum = p->A[0];
+  crossings = p->crossings;
+  {
+    real lat, lon, s12, S12;
+    geod_gendirect(g, p->lat, p->lon, azi, FALSE, s,
+                   &lat, &lon, 0,
+                   0, 0, 0, 0, &S12);
+    tempsum += S12;
+    crossings += transit(p->lon, lon);
+    geod_geninverse(g, lat,  lon, p->lat0,  p->lon0,
+                    &s12, 0, 0, 0, 0, 0, &S12);
+    perimeter += s12;
+    tempsum += S12;
+    crossings += transit(lon, p->lon0);
+  }
+
+  area0 = 4 * pi * g->c2;
+  if (crossings & 1)
+    tempsum += (tempsum < 0 ? 1 : -1) * area0/2;
+  /* area is with the clockwise sense.  If !reverse convert to
+   * counter-clockwise convention. */
+  if (!reverse)
+    tempsum *= -1;
+  /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
+  if (sign) {
+    if (tempsum > area0/2)
+      tempsum -= area0;
+    else if (tempsum <= -area0/2)
+      tempsum += area0;
+  } else {
+    if (tempsum >= area0)
+      tempsum -= area0;
+    else if (tempsum < 0)
+      tempsum += area0;
+  }
+  if (pP) *pP = perimeter;
+  if (pA) *pA = 0 + tempsum;
+  return num;
+}
+
+void geod_polygonarea(const struct geod_geodesic* g,
+                      real lats[], real lons[], int n,
+                      real* pA, real* pP) {
+  int i;
+  struct geod_polygon p;
+  geod_polygon_init(&p, FALSE);
+  for (i = 0; i < n; ++i)
+    geod_polygon_addpoint(g, &p, lats[i], lons[i]);
+  geod_polygon_compute(g, &p, FALSE, TRUE, pA, pP);
+}
+
+/** @endcond */
diff --git a/src/geodesic.h b/src/geodesic.h
index 87d087a5..69985009 100644
--- a/src/geodesic.h
+++ b/src/geodesic.h
@@ -27,8 +27,8 @@
  * - 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
+ * - the inverse problem -- given \e lat1, \e lon1, and \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 \e a (typically in
@@ -97,31 +97,44 @@
  *   azi2] + [\e d, \e d], for arbitrary \e d.
  *
  * 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.
+ * in <a href="http://geographiclib.sf.net"> GeographicLib</a>.  The "class
+ * data" is represented by the structs geod_geodesic, geod_geodesicline,
+ * geod_polygon and pointers to these objects are passed as initial arguments
+ * 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
+ * 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
  * 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.31.
+ * <a href="../index.html">GeographicLib</a> 1.32.
  **********************************************************************/
 
 #if !defined(GEODESIC_H)
 #define GEODESIC_H 1
 
+/**
+ * The major version of the geodesic library.  (This tracks the version of
+ * GeographicLib.)
+ **********************************************************************/
+#define GEODESIC_VERSION_MAJOR 1
+/**
+ * The minor version of the geodesic library.  (This tracks the version of
+ * GeographicLib.)
+ **********************************************************************/
+#define GEODESIC_VERSION_MINOR 32
+/**
+ * The patch level of the geodesic library.  (This tracks the version of
+ * GeographicLib.)
+ **********************************************************************/
+#define GEODESIC_VERSION_PATCH 0
+
 #if defined(__cplusplus)
 extern "C" {
 #endif
@@ -159,30 +172,50 @@ extern "C" {
   };
 
   /**
-   * Initialize a geod_geodesic object
+   * The struct for accumulating information about a geodesic polygon.  This is
+   * used for computing the perimeter and area of a polygon.  This must be
+   * initialized by geod_polygon_init() before use.
+   **********************************************************************/
+  struct geod_polygon {
+    double lat;                 /**< the current latitude */
+    double lon;                 /**< the current longitude */
+    /**< @cond SKIP */
+    double lat0;
+    double lon0;
+    double A[2];
+    double P[2];
+    int polyline;
+    int crossings;
+    /**< @endcond */
+    unsigned num;               /**< the number of points so far */
+  };
+
+  /**
+   * Initialize a geod_geodesic object.
    *
-   * @param[out] g the object to be initialized.
+   * @param[out] g a pointer to 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
+   * Initialize a geod_geodesicline object.
    *
-   * @param[out] l the object to be initialized.
-   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[out] l a pointer to the object to be initialized.
+   * @param[in] g a pointer to 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
+   * @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;).
+   * 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 geod_mask values are [see geod_mask()]:
    * - \e caps |= GEOD_LATITUDE for the latitude \e lat2; this is
@@ -210,7 +243,8 @@ extern "C" {
   /**
    * Solve the direct geodesic problem.
    *
-   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] g a pointer to 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).
@@ -224,7 +258,7 @@ extern "C" {
    * 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
+   * the "return" arguments \e 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
@@ -250,7 +284,8 @@ extern "C" {
   /**
    * Solve the inverse geodesic problem.
    *
-   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] g a pointer to 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).
@@ -264,8 +299,8 @@ extern "C" {
    * 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.
+   * Any of the "return" arguments \e 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, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;), and
@@ -292,18 +327,19 @@ extern "C" {
   /**
    * Compute the position along a geod_geodesicline.
    *
-   * @param[in] l the geod_geodesicline object specifying the geodesic line.
+   * @param[in] l a pointer to 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.
+   *   that \e l 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
+   * \e 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
@@ -340,7 +376,8 @@ extern "C" {
   /**
    * The general direct geodesic problem.
    *
-   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] g a pointer to 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).
@@ -367,7 +404,7 @@ extern "C" {
    * 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
+   * arguments \e plat2, etc., may be replaced by 0, if you do not need some
    * quantities computed.
    **********************************************************************/
   double geod_gendirect(const struct geod_geodesic* g,
@@ -380,7 +417,8 @@ extern "C" {
   /**
    * The general inverse geodesic calculation.
    *
-   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] g a pointer to 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).
@@ -401,7 +439,7 @@ extern "C" {
    * \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
+   * "return" arguments \e ps12, etc., may be replaced by 0, if you do not need
    * some quantities computed.
    **********************************************************************/
   double geod_geninverse(const struct geod_geodesic* g,
@@ -413,7 +451,8 @@ extern "C" {
   /**
    * The general position function.
    *
-   * @param[in] l the geod_geodesicline object specifying the geodesic line.
+   * @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 |=
    *   GEOD_DISTANCE_IN.
@@ -443,11 +482,11 @@ extern "C" {
    * \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
+   * \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.
    *
-   * Example, computer way points between JFK and Singapore Changi Airport
+   * Example, compute 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
@@ -476,9 +515,181 @@ extern "C" {
                           double* pS12);
 
   /**
-   * The area of a geodesic polygon.
+   * Initialize a geod_polygon object.
+   *
+   * @param[out] p a pointer to the object to be initialized.
+   * @param[in] polylinep non-zero if a polyline instead of a polygon.
+   *
+   * If \e polylinep is zero, then the sequence of vertices and edges added by
+   * geod_polygon_addpoint() and geod_polygon_addedge() define a polygon and
+   * the perimeter and area are returned by geod_polygon_compute().  If \e
+   * polylinep is non-zero, then the vertices and edges define a polyline and
+   * only the perimeter is returned by geod_polygon_compute().
+   *
+   * An example of the use of this function is given in the documentation for
+   * geod_polygon_compute().
+   **********************************************************************/
+  void geod_polygon_init(struct geod_polygon* p, int polylinep);
+
+  /**
+   * Add a point to the polygon or polyline.
+   *
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
+   * @param[in,out] p a pointer to the geod_polygon object specifying the
+   *   polygon.
+   * @param[in] lat the latitude of the point (degrees).
+   * @param[in] lon the longitude of the point (degrees).
+   *
+   * \e g and \e p must have been initialized with calls to geod_init() and
+   * geod_polygon_init(), respectively.  The same \e g must be used for all the
+   * points and edges in a polygon.  \e lat should be in the range
+   * [&minus;90&deg;, 90&deg;] and \e lon should be in the range
+   * [&minus;540&deg;, 540&deg;).
+   *
+   * An example of the use of this function is given in the documentation for
+   * geod_polygon_compute().
+   **********************************************************************/
+  void geod_polygon_addpoint(const struct geod_geodesic* g,
+                             struct geod_polygon* p,
+                             double lat, double lon);
+
+  /**
+   * Add an edge to the polygon or polyline.
+   *
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
+   * @param[in,out] p a pointer to the geod_polygon object specifying the
+   *   polygon.
+   * @param[in] azi azimuth at current point (degrees).
+   * @param[in] s distance from current point to next point (meters).
+   *
+   * \e g and \e p must have been initialized with calls to geod_init() and
+   * geod_polygon_init(), respectively.  The same \e g must be used for all the
+   * points and edges in a polygon.  \e azi should be in the range
+   * [&minus;540&deg;, 540&deg;).  This does nothing if no points have been
+   * added yet.  The \e lat and \e lon fields of \e p give the location of
+   * the new vertex.
+   **********************************************************************/
+  void geod_polygon_addedge(const struct geod_geodesic* g,
+                            struct geod_polygon* p,
+                            double azi, double s);
+
+  /**
+   * Return the results for a polygon.
+   *
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
+   * @param[in] p a pointer to the geod_polygon object specifying the polygon.
+   * @param[in] reverse if non-zero then clockwise (instead of
+   *   counter-clockwise) traversal counts as a positive area.
+   * @param[in] sign if non-zero then return a signed result for the area if
+   *   the polygon is traversed in the "wrong" direction instead of returning
+   *   the area for the rest of the earth.
+   * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
+   *   only set if \e polyline is non-zero in the call to geod_polygon_init().
+   * @param[out] pP pointer to the perimeter of the polygon or length of the
+   *   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.
+   *
+   * 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).
+   @code
+   double A, P;
+   int n;
+   struct geod_geodesic g;
+   struct geod_polygon p;
+   geod_init(&g, 6378137, 1/298.257223563);
+   geod_polygon_init(&p, 0);
+
+   geod_polygon_addpoint(&g, &p,  0,  0);
+   geod_polygon_addpoint(&g, &p,  0, 90);
+   geod_polygon_addpoint(&g, &p, 90,  0);
+   n = geod_polygon_compute(&g, &p, 0, 1, &A, &P);
+   printf("%d %.8f %.3f\n", n, P, A);
+   @endcode
+   **********************************************************************/
+  unsigned geod_polygon_compute(const struct geod_geodesic* g,
+                                const struct geod_polygon* p,
+                                int reverse, int sign,
+                                double* pA, double* pP);
+
+  /**
+   * Return the results assuming a tentative final test point is added;
+   * however, the data for the test point is not saved.  This lets you report a
+   * running result for the perimeter and area as the user moves the mouse
+   * cursor.  Ordinary floating point arithmetic is used to accumulate the data
+   * for the test point; thus the area and perimeter returned are less accurate
+   * than if geod_polygon_addpoint() and geod_polygon_compute() are used.
+   *
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
+   * @param[in] p a pointer to the geod_polygon object specifying the polygon.
+   * @param[in] lat the latitude of the test point (degrees).
+   * @param[in] lon the longitude of the test point (degrees).
+   * @param[in] reverse if non-zero then clockwise (instead of
+   *   counter-clockwise) traversal counts as a positive area.
+   * @param[in] sign if non-zero then return a signed result for the area if
+   *   the polygon is traversed in the "wrong" direction instead of returning
+   *   the area for the rest of the earth.
+   * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
+   *   only set if \e polyline is non-zero in the call to geod_polygon_init().
+   * @param[out] pP pointer to the perimeter of the polygon or length of the
+   *   polyline (meters).
+   * @return the number of points.
+   *
+   * \e lat should be in the range [&minus;90&deg;, 90&deg;] and \e
+   * lon should be in the range [&minus;540&deg;, 540&deg;).
+   **********************************************************************/
+  unsigned geod_polygon_testpoint(const struct geod_geodesic* g,
+                                  const struct geod_polygon* p,
+                                  double lat, double lon,
+                                  int reverse, int sign,
+                                  double* pA, double* pP);
+
+  /**
+   * Return the results assuming a tentative final test point is added via an
+   * azimuth and distance; however, the data for the test point is not saved.
+   * This lets you report a running result for the perimeter and area as the
+   * user moves the mouse cursor.  Ordinary floating point arithmetic is used
+   * to accumulate the data for the test point; thus the area and perimeter
+   * returned are less accurate than if geod_polygon_addedge() and
+   * geod_polygon_compute() are used.
+   *
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
+   * @param[in] p a pointer to the geod_polygon object specifying the polygon.
+   * @param[in] azi azimuth at current point (degrees).
+   * @param[in] s distance from current point to final test point (meters).
+   * @param[in] reverse if non-zero then clockwise (instead of
+   *   counter-clockwise) traversal counts as a positive area.
+   * @param[in] sign if non-zero then return a signed result for the area if
+   *   the polygon is traversed in the "wrong" direction instead of returning
+   *   the area for the rest of the earth.
+   * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
+   *   only set if \e polyline is non-zero in the call to geod_polygon_init().
+   * @param[out] pP pointer to the perimeter of the polygon or length of the
+   *   polyline (meters).
+   * @return the number of points.
+   *
+   * \e azi should be in the range [&minus;540&deg;, 540&deg;).
+   **********************************************************************/
+  unsigned geod_polygon_testedge(const struct geod_geodesic* g,
+                                 const struct geod_polygon* p,
+                                 double azi, double s,
+                                 int reverse, int sign,
+                                 double* pA, double* pP);
+
+  /**
+   * A simple interface for computing the area of a geodesic polygon.
    *
-   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] g a pointer to 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.
@@ -512,7 +723,7 @@ extern "C" {
                         double* pA, double* pP);
 
   /**
-   * mask values for geod_geodesicline capabilities.
+   * mask values for the the \e caps argument to geod_lineinit().
    **********************************************************************/
   enum geod_mask {
     GEOD_NONE         = 0U,                     /**< Calculate nothing */