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

12 files changed, 799 insertions, 397 deletions

diff --git a/HOWTO-RELEASE b/HOWTO-RELEASE
index 8ae52131..41bc3090 100644
--- a/HOWTO-RELEASE
+++ b/HOWTO-RELEASE
@@ -22,7 +22,7 @@
    - If any interfaces have been removed since the last public release, then 
      set age to 0. 
 
-4.5) Run "autgen.sh" (hopefully on the same machine it was last run on)
+4.5) Run "autogen.sh" (hopefully on the same machine it was last run on)
 
 5) Add a note to the ChangeLog that a new release is being issued, and what
    the release number is.
diff --git a/man/man1/geod.1 b/man/man1/geod.1
index ab6c74f7..0a2522ea 100644
--- a/man/man1/geod.1
+++ b/man/man1/geod.1
@@ -1,8 +1,8 @@
-.\" @(#)geod.1 - 1.1
-.nr LL 5.5i
+.\" @(#)geod.1
+.nr LL 6.5i
 .ad b
 .hy 1
-.TH GEOD 1 "2012/11/29 Rel. 4.8"
+.TH GEOD 1 "2013/02/27 Rel. 4.8"
 .SH NAME
 geod \- direct geodesic computations
 .br
@@ -171,7 +171,7 @@ and number of points to be determined.
 .RE
 .SH EXAMPLE
 The following script determines the geodesic azimuths and distance in
-U.S. stature miles from Boston, MA, to Portland, OR:
+U.S. statute miles from Boston, MA, to Portland, OR:
 .RS 5
  \f(CWgeod +ellps=clrk66 <<EOF -I +units=us-mi
  42d15'N 71d07'W 45d31'N 123d41'W
@@ -199,9 +199,19 @@ which gives:
 Note: lack of precision in the distance value compromises
 the precision of the Portland location.
 .SH SEE ALSO
-C. F. F. Karney,
-.I "Algorithms for Geodesics," J. Geodesy (2012);
-DOI: 10.1007/s00190-012-0578-z,
+.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
+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,
+.br
 http://geographiclib.sf.net/geod-addenda.html.
 .SH HOME PAGE
-http://www.remotesensing.org/proj
+http://proj.osgeo.org
diff --git a/man/man3/Makefile.am b/man/man3/Makefile.am
index 2a41471a..4a3128a8 100644
--- a/man/man3/Makefile.am
+++ b/man/man3/Makefile.am
@@ -1,3 +1,3 @@
-man_MANS = 	pj_init.3
+man_MANS = 	pj_init.3 geodesic.3
 
 EXTRA_DIST	=	$(man_MANS)
diff --git a/man/man3/Makefile.in b/man/man3/Makefile.in
index faa81dea..c9e3da99 100644
--- a/man/man3/Makefile.in
+++ b/man/man3/Makefile.in
@@ -187,7 +187,7 @@ target_alias = @target_alias@
 top_build_prefix = @top_build_prefix@
 top_builddir = @top_builddir@
 top_srcdir = @top_srcdir@
-man_MANS = pj_init.3
+man_MANS = pj_init.3 geodesic.3
 EXTRA_DIST = $(man_MANS)
 all: all-am
 
diff --git a/src/Makefile.am b/src/Makefile.am
index bbf28ba6..c1d1998a 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -4,15 +4,14 @@ EXTRA_PROGRAMS = multistresstest
 INCLUDES =	-DPROJ_LIB=\"$(pkgdatadir)\" \
 		-DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@
 
-include_HEADERS = proj_api.h org_proj4_Projections.h org_proj4_PJ.h
+include_HEADERS = proj_api.h org_proj4_Projections.h org_proj4_PJ.h geodesic.h
 
 EXTRA_DIST = makefile.vc proj.def
 
 proj_SOURCES = proj.c gen_cheb.c p_series.c
 cs2cs_SOURCES = cs2cs.c gen_cheb.c p_series.c
 nad2bin_SOURCES = nad2bin.c
-geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h \
-	geodesic.c geodesic.h
+geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h
 multistresstest_SOURCES = multistresstest.c
 
 proj_LDADD = libproj.la
@@ -64,8 +63,7 @@ libproj_la_SOURCES = \
 	nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \
 	pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \
 	geocent.c geocent.h pj_utils.c pj_gridinfo.c pj_gridlist.c \
-	jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c
-
+	jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c geodesic.c
 
 install-exec-local:
 	rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT)
diff --git a/src/Makefile.in b/src/Makefile.in
index eae1c90c..dd32d94b 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -107,7 +107,7 @@ am_libproj_la_OBJECTS = PJ_aeqd.lo PJ_gnom.lo PJ_laea.lo \
 	pj_apply_gridshift.lo pj_datums.lo pj_datum_set.lo \
 	pj_transform.lo geocent.lo pj_utils.lo pj_gridinfo.lo \
 	pj_gridlist.lo jniproj.lo pj_mutex.lo pj_initcache.lo \
-	pj_apply_vgridshift.lo
+	pj_apply_vgridshift.lo geodesic.lo
 libproj_la_OBJECTS = $(am_libproj_la_OBJECTS)
 libproj_la_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
 	$(LIBTOOLFLAGS) --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \
@@ -118,7 +118,7 @@ am_cs2cs_OBJECTS = cs2cs.$(OBJEXT) gen_cheb.$(OBJEXT) \
 cs2cs_OBJECTS = $(am_cs2cs_OBJECTS)
 cs2cs_DEPENDENCIES = libproj.la
 am_geod_OBJECTS = geod.$(OBJEXT) geod_set.$(OBJEXT) \
-	geod_interface.$(OBJEXT) geodesic.$(OBJEXT)
+	geod_interface.$(OBJEXT)
 geod_OBJECTS = $(am_geod_OBJECTS)
 geod_DEPENDENCIES = libproj.la
 am_multistresstest_OBJECTS = multistresstest.$(OBJEXT)
@@ -270,14 +270,12 @@ top_srcdir = @top_srcdir@
 INCLUDES = -DPROJ_LIB=\"$(pkgdatadir)\" \
 		-DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@
 
-include_HEADERS = proj_api.h org_proj4_Projections.h org_proj4_PJ.h
+include_HEADERS = proj_api.h org_proj4_Projections.h org_proj4_PJ.h geodesic.h
 EXTRA_DIST = makefile.vc proj.def
 proj_SOURCES = proj.c gen_cheb.c p_series.c
 cs2cs_SOURCES = cs2cs.c gen_cheb.c p_series.c
 nad2bin_SOURCES = nad2bin.c
-geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h \
-	geodesic.c geodesic.h
-
+geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h
 multistresstest_SOURCES = multistresstest.c
 proj_LDADD = libproj.la
 cs2cs_LDADD = libproj.la
@@ -325,7 +323,7 @@ libproj_la_SOURCES = \
 	nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \
 	pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \
 	geocent.c geocent.h pj_utils.c pj_gridinfo.c pj_gridlist.c \
-	jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c
+	jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c geodesic.c
 
 all: proj_config.h
 	$(MAKE) $(AM_MAKEFLAGS) all-am
@@ -583,7 +581,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geod.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geod_interface.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geod_set.Po@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geodesic.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geodesic.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/jniproj.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/mk_cheby.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/multistresstest.Po@am__quote@
diff --git a/src/geod_interface.c b/src/geod_interface.c
index c144a35d..dc757282 100644
--- a/src/geod_interface.c
+++ b/src/geod_interface.c
@@ -2,20 +2,20 @@
 #include "geod_interface.h"
 
 void geod_ini(void) {
-  GeodesicInit(&GlobalGeodesic, geod_a, geod_f);
+  geod_init(&GlobalGeodesic, geod_a, geod_f);
 }
 
 void geod_pre(void) {
   double
     degree = PI/180,
     lat1 = phi1 / degree, lon1 = lam1 /degree, azi1 = al12 / degree;
-  GeodesicLineInit(&GlobalGeodesicLine, &GlobalGeodesic,
+  geod_lineinit(&GlobalGeodesicLine, &GlobalGeodesic,
 		   lat1, lon1, azi1, 0U);
 }
 
 void geod_for(void) {
   double degree = PI/180, s12 = geod_S, lat2, lon2, azi2;
-  Position(&GlobalGeodesicLine, s12, &lat2, &lon2, &azi2);
+  geod_position(&GlobalGeodesicLine, s12, &lat2, &lon2, &azi2);
   azi2 += azi2 >= 0 ? -180 : 180; /* Compute back azimuth */
   phi2 = lat2 * degree;
   lam2 = lon2 * degree;
@@ -28,7 +28,7 @@ void geod_inv(void) {
     lat1 = phi1 / degree, lon1 = lam1 / degree,
     lat2 = phi2 / degree, lon2 = lam2 / degree,
     azi1, azi2, s12;
-  Inverse(&GlobalGeodesic, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2);
+  geod_inverse(&GlobalGeodesic, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2);
   azi2 += azi2 >= 0 ? -180 : 180; /* Compute back azimuth */
   al12 = azi1 * degree; al21 = azi2 * degree; geod_S = s12;
 }
diff --git a/src/geod_interface.h b/src/geod_interface.h
index 161c7572..255d505a 100644
--- a/src/geod_interface.h
+++ b/src/geod_interface.h
@@ -27,8 +27,8 @@ GEOD_EXTERN struct geodesic {
 # define al21	GEODESIC.ALPHA21
 # define geod_S	GEODESIC.DIST
     
-GEOD_EXTERN struct Geodesic GlobalGeodesic;
-GEOD_EXTERN struct GeodesicLine GlobalGeodesicLine;
+GEOD_EXTERN struct geod_geodesic GlobalGeodesic;
+GEOD_EXTERN struct geod_geodesicline GlobalGeodesicLine;
 GEOD_EXTERN int n_alpha, n_S;
 GEOD_EXTERN double to_meter, fr_meter, del_alpha;
 	
diff --git a/src/geodesic.c b/src/geodesic.c
index 794439a6..c5d481ee 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -1,22 +1,27 @@
+/**
+ * \file geodesic.c
+ * \brief Implementation of the geodesic routines in C
+ **********************************************************************/
+
 /*
  * This is a C implementation of the geodesic algorithms described in
  *
- *   C. F. F. Karney
- *   Algorithms for geodesics
- *   J. Geodesy (2012)
+ *   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
  *   Addenda: http://geographiclib.sf.net/geod-addenda.html
  *
  * See the comments in geodesic.h for documentation.
  *
- * 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.
  */
 
 #include "geodesic.h"
+
+/** @cond SKIP */
 #include <math.h>
 
 #define GEOGRAPHICLIB_GEODESIC_ORDER 6
@@ -117,6 +122,19 @@ static real cbrtx(real x) {
   return x < 0 ? -y : y;
 }
 
+static real sumx(real u, real v, real* t) {
+  volatile real s = u + v;
+  volatile real up = s - v;
+  volatile real vpp = s - up;
+  up -= u;
+  vpp -= v;
+  *t = -(up + vpp);
+  /* error-free sum:
+   * u + v =       s      + t
+   *       = round(u + v) + t */
+  return s;
+}
+
 static real minx(real x, real y)
 { return x < y ? x : y; }
 
@@ -137,21 +155,30 @@ static real AngNormalize(real x)
 static real AngNormalize2(real x)
 { return AngNormalize(fmod(x, (real)(360))); }
 
+static real AngDiff(real x, real y) {
+  real t, d = sumx(-x, y, &t);
+  if ((d - (real)(180)) + t > (real)(0))       /* y - x > 180 */
+    d -= (real)(360);                          /* exact */
+  else if ((d + (real)(180)) + t <= (real)(0)) /* y - x <= -180 */
+    d += (real)(360);                          /* exact */
+  return d + t;
+}
+
 static real AngRound(real x) {
-  const real z = (real)(0.0625); /* 1/16 */
+  const real z = 1/(real)(16);
   volatile real y = fabs(x);
   /* The compiler mustn't "simplify" z - (z - y) to y */
   y = y < z ? z - (z - y) : y;
   return x < 0 ? -y : y;
 }
 
-static void A3coeff(struct Geodesic* g);
-static void C3coeff(struct Geodesic* g);
-static void C4coeff(struct Geodesic* g);
+static void A3coeff(struct geod_geodesic* g);
+static void C3coeff(struct geod_geodesic* g);
+static void C4coeff(struct geod_geodesic* g);
 static real SinCosSeries(boolx sinp,
                          real sinx, real cosx,
                          const real c[], int n);
-static void Lengths(const struct Geodesic* g,
+static void Lengths(const struct geod_geodesic* g,
                     real eps, real sig12,
                     real ssig1, real csig1, real dn1,
                     real ssig2, real csig2, real dn2,
@@ -161,7 +188,7 @@ static void Lengths(const struct Geodesic* g,
                     /* Scratch areas of the right size */
                     real C1a[], real C2a[]);
 static real Astroid(real x, real y);
-static real InverseStart(const struct Geodesic* g,
+static real InverseStart(const struct geod_geodesic* g,
                          real sbet1, real cbet1, real dn1,
                          real sbet2, real cbet2, real dn2,
                          real lam12,
@@ -170,7 +197,7 @@ static real InverseStart(const struct Geodesic* g,
                          real* psalp2, real* pcalp2,
                          /* Scratch areas of the right size */
                          real C1a[], real C2a[]);
-static real Lambda12(const struct Geodesic* g,
+static real Lambda12(const struct geod_geodesic* g,
                      real sbet1, real cbet1, real dn1,
                      real sbet2, real cbet2, real dn2,
                      real salp1, real calp1,
@@ -182,16 +209,19 @@ static real Lambda12(const struct Geodesic* g,
                      boolx diffp, real* pdlam12,
                      /* Scratch areas of the right size */
                      real C1a[], real C2a[], real C3a[]);
-static real A3f(const struct Geodesic* g, real eps);
-static void C3f(const struct Geodesic* g, real eps, real c[]);
-static void C4f(const struct Geodesic* g, real eps, real c[]);
+static real A3f(const struct geod_geodesic* g, real eps);
+static void C3f(const struct geod_geodesic* g, real eps, real c[]);
+static void C4f(const struct geod_geodesic* g, real eps, real c[]);
 static real A1m1f(real eps);
 static void C1f(real eps, real c[]);
 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);
+
+/** @endcond */
 
-void GeodesicInit(struct Geodesic* g, real a, real f) {
+void geod_init(struct geod_geodesic* g, real a, real f) {
   if (!init) Init();
   g->a = a;
   g->f = f <= 1 ? f : 1/f;
@@ -211,9 +241,9 @@ void GeodesicInit(struct Geodesic* g, real a, real f) {
   C4coeff(g);
 }
 
-void GeodesicLineInit(struct GeodesicLine* l,
-                      const struct Geodesic* g,
-                      real lat1, real lon1, real azi1, unsigned caps) {
+void geod_lineinit(struct geod_geodesicline* l,
+                   const struct geod_geodesic* g,
+                   real lat1, real lon1, real azi1, unsigned caps) {
   real alp1, cbet1, sbet1, phi, eps;
   l->a = g->a;
   l->f = g->f;
@@ -221,23 +251,12 @@ void GeodesicLineInit(struct GeodesicLine* l,
   l->c2 = g->c2;
   l->f1 = g->f1;
   /* If caps is 0 assume the standard direct calculation */
-  l->caps = (caps ? caps : DISTANCE_IN | LONGITUDE) |
-    LATITUDE | AZIMUTH; /* Always allow latitude and azimuth */
+  l->caps = (caps ? caps : GEOD_DISTANCE_IN | GEOD_LONGITUDE) |
+    GEOD_LATITUDE | GEOD_AZIMUTH; /* Always allow latitude and azimuth */
 
-  azi1 = AngNormalize(azi1);
+  /* Guard against underflow in salp0 */
+  azi1 = AngRound(AngNormalize(azi1));
   lon1 = AngNormalize(lon1);
-  if (lat1 == 90) {
-    lon1 += lon1 < 0 ? 180 : -180;
-    lon1 = AngNormalize(lon1 - azi1);
-    azi1 = -180;
-  } else if (lat1 == -90) {
-    lon1 = AngNormalize(lon1 + azi1);
-    azi1 = 0;
-  } else {
-    /* Guard against underflow in salp0 */
-    azi1 = AngRound(azi1);
-  }
-
   l->lat1 = lat1;
   l->lon1 = lon1;
   l->azi1 = azi1;
@@ -245,7 +264,7 @@ void GeodesicLineInit(struct GeodesicLine* l,
   alp1 = 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->salp1 =      azi1  == -180 ? 0 : sin(alp1);
   l->calp1 = fabs(azi1) ==   90 ? 0 : cos(alp1);
   phi = lat1 * degree;
   /* Ensure cbet1 = +epsilon at poles */
@@ -312,12 +331,12 @@ void GeodesicLineInit(struct GeodesicLine* l,
   }
 }
 
-real GenPosition(const struct GeodesicLine* l,
-                 boolx arcmode, real s12_a12,
-                 real* plat2, real* plon2, real* pazi2,
-                 real* ps12, real* pm12,
-                 real* pM12, real* pM21,
-                 real* pS12) {
+real geod_genposition(const struct geod_geodesicline* l,
+                      boolx arcmode, real s12_a12,
+                      real* plat2, real* plon2, real* pazi2,
+                      real* ps12, real* pm12,
+                      real* pM12, real* pM21,
+                      real* pS12) {
   real lat2 = 0, lon2 = 0, azi2 = 0, s12 = 0,
     m12 = 0, M12 = 0, M21 = 0, S12 = 0;
   /* Avoid warning about uninitialized B12. */
@@ -325,16 +344,17 @@ real GenPosition(const struct GeodesicLine* l,
   real omg12, lam12, lon12;
   real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2, dn2;
   unsigned outmask =
-    (plat2 ? LATITUDE : 0U) |
-    (plon2 ? LONGITUDE : 0U) |
-    (pazi2 ? AZIMUTH : 0U) |
-    (ps12 ? DISTANCE : 0U) |
-    (pm12 ? REDUCEDLENGTH : 0U) |
-    (pM12 || pM21 ? GEODESICSCALE : 0U) |
-    (pS12 ? AREA : 0U);
+    (plat2 ? GEOD_LATITUDE : 0U) |
+    (plon2 ? GEOD_LONGITUDE : 0U) |
+    (pazi2 ? GEOD_AZIMUTH : 0U) |
+    (ps12 ? GEOD_DISTANCE : 0U) |
+    (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
+    (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
+    (pS12 ? GEOD_AREA : 0U);
 
   outmask &= l->caps & OUT_ALL;
-  if (!( TRUE /*Init()*/ && (arcmode || (l->caps & DISTANCE_IN & OUT_ALL)) ))
+  if (!( TRUE /*Init()*/ &&
+         (arcmode || (l->caps & GEOD_DISTANCE_IN & OUT_ALL)) ))
     /* Uninitialized or impossible distance calculation requested */
     return NaN;
 
@@ -397,7 +417,7 @@ real GenPosition(const struct GeodesicLine* l,
   ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12;
   csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;
   dn2 = sqrt(1 + l->k2 * sq(ssig2));
-  if (outmask & (DISTANCE | REDUCEDLENGTH | GEODESICSCALE)) {
+  if (outmask & (GEOD_DISTANCE | GEOD_REDUCEDLENGTH | GEOD_GEODESICSCALE)) {
     if (arcmode || fabs(l->f) > 0.01)
       B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
     AB1 = (1 + l->A1m1) * (B12 - l->B11);
@@ -417,10 +437,10 @@ real GenPosition(const struct GeodesicLine* l,
   omg12 = atan2(somg2 * l->comg1 - comg2 * l->somg1,
                 comg2 * l->comg1 + somg2 * l->somg1);
 
-  if (outmask & DISTANCE)
+  if (outmask & GEOD_DISTANCE)
     s12 = arcmode ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12;
 
-  if (outmask & LONGITUDE) {
+  if (outmask & GEOD_LONGITUDE) {
     lam12 = omg12 + l->A3c *
       ( sig12 + (SinCosSeries(TRUE, ssig2, csig2, l->C3a, nC3-1)
                  - l->B31));
@@ -431,31 +451,31 @@ real GenPosition(const struct GeodesicLine* l,
     lon2 = AngNormalize(l->lon1 + lon12);
   }
 
-  if (outmask & LATITUDE)
+  if (outmask & GEOD_LATITUDE)
     lat2 = atan2(sbet2, l->f1 * cbet2) / degree;
 
-  if (outmask & AZIMUTH)
+  if (outmask & GEOD_AZIMUTH)
     /* minus signs give range [-180, 180). 0- converts -0 to +0. */
     azi2 = 0 - atan2(-salp2, calp2) / degree;
 
-  if (outmask & (REDUCEDLENGTH | GEODESICSCALE)) {
+  if (outmask & (GEOD_REDUCEDLENGTH | GEOD_GEODESICSCALE)) {
     real
       B22 = SinCosSeries(TRUE, ssig2, csig2, l->C2a, nC2),
       AB2 = (1 + l->A2m1) * (B22 - l->B21),
       J12 = (l->A1m1 - l->A2m1) * sig12 + (AB1 - AB2);
-    if (outmask & REDUCEDLENGTH)
+    if (outmask & GEOD_REDUCEDLENGTH)
       /* Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
        * accurate cancellation in the case of coincident points. */
       m12 = l->b * ((dn2 * (l->csig1 * ssig2) - l->dn1 * (l->ssig1 * csig2))
                     - l->csig1 * csig2 * J12);
-    if (outmask & GEODESICSCALE) {
+    if (outmask & GEOD_GEODESICSCALE) {
       real t = l->k2 * (ssig2 - l->ssig1) * (ssig2 + l->ssig1) / (l->dn1 + dn2);
       M12 = csig12 + (t *  ssig2 -  csig2 * J12) * l->ssig1 / l->dn1;
       M21 = csig12 - (t * l->ssig1 - l->csig1 * J12) *  ssig2 /  dn2;
     }
   }
 
-  if (outmask & AREA) {
+  if (outmask & GEOD_AREA) {
     real
       B42 = SinCosSeries(FALSE, ssig2, csig2, l->C4a, nC4);
     real salp12, calp12;
@@ -488,66 +508,66 @@ real GenPosition(const struct GeodesicLine* l,
     S12 = l->c2 * atan2(salp12, calp12) + l->A4 * (B42 - l->B41);
   }
 
-  if (outmask & LATITUDE)
+  if (outmask & GEOD_LATITUDE)
     *plat2 = lat2;
-  if (outmask & LONGITUDE)
+  if (outmask & GEOD_LONGITUDE)
     *plon2 = lon2;
-  if (outmask & AZIMUTH)
+  if (outmask & GEOD_AZIMUTH)
     *pazi2 = azi2;
-  if (outmask & DISTANCE)
+  if (outmask & GEOD_DISTANCE)
     *ps12 = s12;
-  if (outmask & REDUCEDLENGTH)
+  if (outmask & GEOD_REDUCEDLENGTH)
     *pm12 = m12;
-  if (outmask & GEODESICSCALE) {
+  if (outmask & GEOD_GEODESICSCALE) {
     if (pM12) *pM12 = M12;
     if (pM21) *pM21 = M21;
   }
-  if (outmask & AREA)
+  if (outmask & GEOD_AREA)
     *pS12 = S12;
 
   return arcmode ? s12_a12 : sig12 / degree;
 }
 
-void Position(const struct GeodesicLine* l, real s12,
-              real* plat2, real* plon2, real* pazi2) {
-  GenPosition(l, FALSE, s12, plat2, plon2, pazi2, 0, 0, 0, 0, 0);
+void geod_position(const struct geod_geodesicline* l, real s12,
+                   real* plat2, real* plon2, real* pazi2) {
+  geod_genposition(l, FALSE, s12, plat2, plon2, pazi2, 0, 0, 0, 0, 0);
 }
 
-real GenDirect(const struct Geodesic* g,
-               real lat1, real lon1, real azi1,
-               boolx arcmode, real s12_a12,
-               real* plat2, real* plon2, real* pazi2,
-               real* ps12, real* pm12, real* pM12, real* pM21,
-               real* pS12) {
-  struct GeodesicLine l;
+real geod_gendirect(const struct geod_geodesic* g,
+                    real lat1, real lon1, real azi1,
+                    boolx arcmode, real s12_a12,
+                    real* plat2, real* plon2, real* pazi2,
+                    real* ps12, real* pm12, real* pM12, real* pM21,
+                    real* pS12) {
+  struct geod_geodesicline l;
   unsigned outmask =
-    (plat2 ? LATITUDE : 0U) |
-    (plon2 ? LONGITUDE : 0U) |
-    (pazi2 ? AZIMUTH : 0U) |
-    (ps12 ? DISTANCE : 0U) |
-    (pm12 ? REDUCEDLENGTH : 0U) |
-    (pM12 || pM21 ? GEODESICSCALE : 0U) |
-    (pS12 ? AREA : 0U);
-
-  GeodesicLineInit(&l, g, lat1, lon1, azi1,
-                   /* Automatically supply DISTANCE_IN if necessary */
-                   outmask | (arcmode ? NONE : DISTANCE_IN));
-  return GenPosition(&l, arcmode, s12_a12,
-                     plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12);
+    (plat2 ? GEOD_LATITUDE : 0U) |
+    (plon2 ? GEOD_LONGITUDE : 0U) |
+    (pazi2 ? GEOD_AZIMUTH : 0U) |
+    (ps12 ? GEOD_DISTANCE : 0U) |
+    (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
+    (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
+    (pS12 ? GEOD_AREA : 0U);
+
+  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,
+                          plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12);
 }
 
-void Direct(const struct Geodesic* g,
-            real lat1, real lon1, real azi1,
-            real s12,
-            real* plat2, real* plon2, real* pazi2) {
-  GenDirect(g, lat1, lon1, azi1, FALSE, s12, plat2, plon2, pazi2,
-            0, 0, 0, 0, 0);
+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,
+                 0, 0, 0, 0, 0);
 }
 
-real GenInverse(const struct Geodesic* g,
-                real lat1, real lon1, real lat2, real lon2,
-                real* ps12, real* pazi1, real* pazi2,
-                real* pm12, real* pM12, real* pM21, real* pS12) {
+real geod_geninverse(const struct geod_geodesic* g,
+                     real lat1, real lon1, real lat2, real lon2,
+                     real* ps12, real* pazi1, real* pazi2,
+                     real* pm12, real* pM12, real* pM21, real* pS12) {
   real s12 = 0, azi1 = 0, azi2 = 0, m12 = 0, M12 = 0, M21 = 0, S12 = 0;
   real lon12;
   int latsign, lonsign, swapp;
@@ -560,23 +580,22 @@ real GenInverse(const struct Geodesic* g,
   real omg12 = 0;
 
   unsigned outmask =
-    (ps12 ? DISTANCE : 0U) |
-    (pazi1 || pazi2 ? AZIMUTH : 0U) |
-    (pm12 ? REDUCEDLENGTH : 0U) |
-    (pM12 || pM21 ? GEODESICSCALE : 0U) |
-    (pS12 ? AREA : 0U);
+    (ps12 ? GEOD_DISTANCE : 0U) |
+    (pazi1 || pazi2 ? GEOD_AZIMUTH : 0U) |
+    (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
+    (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
+    (pS12 ? GEOD_AREA : 0U);
 
   outmask &= OUT_ALL;
-  lon12 = AngNormalize(AngNormalize(lon2) -
-                            AngNormalize(lon1));
-  /* If very close to being on the same meridian, then make it so.
-   * Not sure this is necessary... */
+  /* Compute longitude difference (AngDiff does this carefully).  Result is
+   * in [-180, 180] but -180 is only for west-going geodesics.  180 is for
+   * east-going and meridional geodesics. */
+  lon12 = AngDiff(AngNormalize(lon1), AngNormalize(lon2));
+  /* If very close to being on the same half-meridian, then make it so. */
   lon12 = AngRound(lon12);
   /* Make longitude difference positive. */
   lonsign = lon12 >= 0 ? 1 : -1;
   lon12 *= lonsign;
-  if (lon12 == 180)
-    lonsign = 1;
   /* If really close to the equator, treat as on equator. */
   lat1 = AngRound(lat1);
   lat2 = AngRound(lat2);
@@ -637,7 +656,6 @@ real GenInverse(const struct Geodesic* g,
   slam12 = lon12 == 180 ? 0 : sin(lam12);
   clam12 = cos(lam12);      /* lon12 == 90 isn't interesting */
 
-
   meridian = lat1 == -90 || slam12 == 0;
 
   if (meridian) {
@@ -660,7 +678,7 @@ real GenInverse(const struct Geodesic* g,
       real dummy;
       Lengths(g, g->n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
               cbet1, cbet2, &s12x, &m12x, &dummy,
-              (outmask & GEODESICSCALE) != 0U, &M12, &M21, C1a, C2a);
+              (outmask & GEOD_GEODESICSCALE) != 0U, &M12, &M21, C1a, C2a);
     }
     /* Add the check for sig12 since zero length geodesics might yield m12 <
      * 0.  Test case was
@@ -688,7 +706,7 @@ real GenInverse(const struct Geodesic* g,
     s12x = g->a * lam12;
     sig12 = omg12 = lam12 / g->f1;
     m12x = g->b * sin(sig12);
-    if (outmask & GEODESICSCALE)
+    if (outmask & GEOD_GEODESICSCALE)
       M12 = M21 = cos(sig12);
     a12 = lon12 / g->f1;
 
@@ -708,7 +726,7 @@ real GenInverse(const struct Geodesic* g,
       real dnm = (dn1 + dn2) / 2;
       s12x = sig12 * g->b * dnm;
       m12x = sq(dnm) * g->b * sin(sig12 / dnm);
-      if (outmask & GEODESICSCALE)
+      if (outmask & GEOD_GEODESICSCALE)
         M12 = M21 = cos(sig12 / dnm);
       a12 = sig12 / degree;
       omg12 = lam12 / (g->f1 * dnm);
@@ -739,13 +757,13 @@ real GenInverse(const struct Geodesic* g,
                         &eps, &omg12, numit < maxit1, &dv, C1a, C2a, C3a)
                - lam12);
         /* 2 * tol0 is approximately 1 ulp for a number in [0, pi]. */
-        if (tripb || fabs(v) < (tripn ? 8 : 2) * tol0) break;
+        /* Reversed test to allow escape with NaNs */
+        if (tripb || !(fabs(v) >= (tripn ? 8 : 2) * tol0)) break;
         /* Update bracketing values */
-        if (v > 0 && (numit > maxit1 || calp1/salp1 > calp1b/salp1b)) {
-          salp1b = salp1; calp1b = calp1;
-        } else if (numit > maxit1 || calp1/salp1 < calp1a/salp1a) {
-          salp1a = salp1; calp1a = calp1;
-        }
+        if (v > 0 && (numit > maxit1 || calp1/salp1 > calp1b/salp1b))
+          { salp1b = salp1; calp1b = calp1; }
+        else if (v < 0 && (numit > maxit1 || calp1/salp1 < calp1a/salp1a))
+          { salp1a = salp1; calp1a = calp1; }
         if (numit < maxit1 && dv > 0) {
           real
             dalp1 = -v/dv;
@@ -782,7 +800,7 @@ real GenInverse(const struct Geodesic* g,
         real dummy;
         Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
                 cbet1, cbet2, &s12x, &m12x, &dummy,
-                (outmask & GEODESICSCALE) != 0U, &M12, &M21, C1a, C2a);
+                (outmask & GEOD_GEODESICSCALE) != 0U, &M12, &M21, C1a, C2a);
       }
       m12x *= g->b;
       s12x *= g->b;
@@ -791,13 +809,13 @@ real GenInverse(const struct Geodesic* g,
     }
   }
 
-  if (outmask & DISTANCE)
+  if (outmask & GEOD_DISTANCE)
     s12 = 0 + s12x;             /* Convert -0 to 0 */
 
-  if (outmask & REDUCEDLENGTH)
+  if (outmask & GEOD_REDUCEDLENGTH)
     m12 = 0 + m12x;             /* Convert -0 to 0 */
 
-  if (outmask & AREA) {
+  if (outmask & GEOD_AREA) {
     real
       /* From Lambda12: sin(alp1) * cos(bet1) = sin(alp0) */
       salp0 = salp1 * cbet1,
@@ -860,44 +878,46 @@ real GenInverse(const struct Geodesic* g,
   if (swapp < 0) {
     swapx(&salp1, &salp2);
     swapx(&calp1, &calp2);
-    if (outmask & GEODESICSCALE)
+    if (outmask & GEOD_GEODESICSCALE)
       swapx(&M12, &M21);
   }
 
   salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
   salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
 
-  if (outmask & AZIMUTH) {
+  if (outmask & GEOD_AZIMUTH) {
     /* minus signs give range [-180, 180). 0- converts -0 to +0. */
     azi1 = 0 - atan2(-salp1, calp1) / degree;
     azi2 = 0 - atan2(-salp2, calp2) / degree;
   }
 
-  if (outmask & DISTANCE)
+  if (outmask & GEOD_DISTANCE)
     *ps12 = s12;
-  if (outmask & AZIMUTH) {
+  if (outmask & GEOD_AZIMUTH) {
     if (pazi1) *pazi1 = azi1;
     if (pazi2) *pazi2 = azi2;
   }
-  if (outmask & REDUCEDLENGTH)
+  if (outmask & GEOD_REDUCEDLENGTH)
     *pm12 = m12;
-  if (outmask & GEODESICSCALE) {
+  if (outmask & GEOD_GEODESICSCALE) {
     if (pM12) *pM12 = M12;
     if (pM21) *pM21 = M21;
   }
-  if (outmask & AREA)
+  if (outmask & GEOD_AREA)
     *pS12 = S12;
 
   /* Returned value in [0, 180] */
   return a12;
 }
 
-void Inverse(const struct Geodesic* g,
-             double lat1, double lon1, double lat2, double lon2,
-             double* ps12, double* pazi1, double* pazi2) {
-  GenInverse(g, lat1, lon1, lat2, lon2, ps12, pazi1, pazi2, 0, 0, 0, 0);
+void geod_inverse(const struct geod_geodesic* g,
+                  double lat1, double lon1, double lat2, double lon2,
+                  double* ps12, double* pazi1, double* 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) {
@@ -922,7 +942,7 @@ real SinCosSeries(boolx sinp,
     : cosx * (y0 - y1);         /* cos(x) * (y0 - y1) */
 }
 
-void Lengths(const struct Geodesic* g,
+void Lengths(const struct geod_geodesic* g,
              real eps, real sig12,
              real ssig1, real csig1, real dn1,
              real ssig2, real csig2, real dn2,
@@ -1019,7 +1039,7 @@ real Astroid(real x, real y) {
   return k;
 }
 
-real InverseStart(const struct Geodesic* g,
+real InverseStart(const struct geod_geodesic* g,
                   real sbet1, real cbet1, real dn1,
                   real sbet2, real cbet2, real dn2,
                   real lam12,
@@ -1187,7 +1207,7 @@ real InverseStart(const struct Geodesic* g,
   return sig12;
 }
 
-real Lambda12(const struct Geodesic* g,
+real Lambda12(const struct geod_geodesic* g,
               real sbet1, real cbet1, real dn1,
               real sbet2, real cbet2, real dn2,
               real salp1, real calp1,
@@ -1286,7 +1306,7 @@ real Lambda12(const struct Geodesic* g,
   return lam12;
 }
 
-real A3f(const struct Geodesic* g, real eps) {
+real A3f(const struct geod_geodesic* g, real eps) {
   /* Evaluate sum(A3x[k] * eps^k, k, 0, nA3x-1) by Horner's method */
   real v = 0;
   int i;
@@ -1295,7 +1315,7 @@ real A3f(const struct Geodesic* g, real eps) {
   return v;
 }
 
-void C3f(const struct Geodesic* g, real eps, real c[]) {
+void C3f(const struct geod_geodesic* g, real eps, real c[]) {
   /* Evaluate C3 coeffs by Horner's method
    * Elements c[1] thru c[nC3 - 1] are set */
   int i, j, k;
@@ -1313,7 +1333,7 @@ void C3f(const struct Geodesic* g, real eps, real c[]) {
   }
 }
 
-void C4f(const struct Geodesic* g, real eps, real c[]) {
+void C4f(const struct geod_geodesic* g, real eps, real c[]) {
   /* Evaluate C4 coeffs by Horner's method
    * Elements c[0] thru c[nC4 - 1] are set */
   int i, j, k;
@@ -1404,7 +1424,7 @@ void C2f(real eps, real c[])  {
 }
 
 /* The scale factor A3 = mean value of (d/dsigma)I3 */
-void A3coeff(struct Geodesic* g) {
+void A3coeff(struct geod_geodesic* g) {
   g->A3x[0] = 1;
   g->A3x[1] = (g->n-1)/2;
   g->A3x[2] = (g->n*(3*g->n-1)-2)/8;
@@ -1414,7 +1434,7 @@ void A3coeff(struct Geodesic* g) {
 }
 
 /* The coefficients C3[l] in the Fourier expansion of B3 */
-void C3coeff(struct Geodesic* g) {
+void C3coeff(struct geod_geodesic* g) {
   g->C3x[0] = (1-g->n)/4;
   g->C3x[1] = (1-g->n*g->n)/8;
   g->C3x[2] = ((3-g->n)*g->n+3)/64;
@@ -1435,7 +1455,7 @@ void C3coeff(struct Geodesic* g) {
 /* Generated by Maxima on 2012-10-19 08:02:34-04:00 */
 
 /* The coefficients C4[l] in the Fourier expansion of I4 */
-void C4coeff(struct Geodesic* g) {
+void C4coeff(struct geod_geodesic* g) {
   g->C4x[0] = (g->n*(g->n*(g->n*(g->n*(100*g->n+208)+572)+3432)-12012)+30030)/
     45045;
   g->C4x[1] = (g->n*(g->n*(g->n*(64*g->n+624)-4576)+6864)-3003)/15015;
@@ -1459,3 +1479,41 @@ void C4coeff(struct Geodesic* g) {
   g->C4x[19] = -128/(real)(135135);
   g->C4x[20] = 128/(real)(99099);
 }
+
+int transit(real lon1, real lon2) {
+  real lon12;
+  /* Return 1 or -1 if crossing prime meridian in east or west direction.
+   * Otherwise return zero. */
+  /* Compute lon12 the same way as Geodesic::Inverse. */
+  lon1 = AngNormalize(lon1);
+  lon2 = AngNormalize(lon2);
+  lon12 = AngDiff(lon1, lon2);
+  return lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 :
+    (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0);
+}
+
+/** @endcond */
+
+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 = 0, P = 0;
+  for (i = 0; i < n; ++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);
+    P += s12;
+    A -= S12; /* The minus sign is due to the counter-clockwise convention */
+    crossings += transit(lons[i], lons[(i + 1) % n]);
+  }
+  if (crossings & 1)
+    A += (A < 0 ? 1 : -1) * area0/2;
+  /* Put area in (-area0/2, area0/2] */
+  if (A > area0/2)
+    A -= area0;
+  else if (A <= -area0/2)
+    A += area0;
+  if (pA) *pA = A;
+  if (pP) *pP = P;
+}
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)
diff --git a/src/makefile.vc b/src/makefile.vc
index c9b68cd5..bd042332 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -55,11 +55,12 @@ support = \
     pj_utils.obj pj_gridlist.obj pj_gridinfo.obj \
 	proj_mdist.obj pj_mutex.obj pj_initcache.obj \
 	pj_ctx.obj pj_log.obj pj_apply_vgridshift.obj
-
-LIBOBJ	=	$(support) $(pseudo) $(azimuthal) $(conic) $(cylinder) $(misc)
+geodesic = geodesic.obj
+LIBOBJ	=	$(support) $(pseudo) $(azimuthal) $(conic) $(cylinder) $(misc) \
+	$(geodesic)
 PROJEXE_OBJ	= proj.obj gen_cheb.obj p_series.obj emess.obj
 CS2CSEXE_OBJ	= cs2cs.obj gen_cheb.obj p_series.obj emess.obj
-GEODEXE_OBJ	= geod.obj geod_set.obj geod_interface.obj geodesic.obj emess.obj
+GEODEXE_OBJ	= geod.obj geod_set.obj geod_interface.obj emess.obj
 PROJ_DLL 	= proj$(VERSION).dll
 PROJ_EXE    = proj.exe
 CS2CS_EXE   = cs2cs.exe
@@ -134,4 +135,5 @@ install: all
 	copy *.dll $(INSTDIR)\bin
 	copy *.lib $(INSTDIR)\lib
 	copy proj_api.h $(INSTDIR)\include
+	copy geodesic.h $(INSTDIR)\include
 	
diff --git a/src/proj.def b/src/proj.def
index f4d60cab..7ae3daaa 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -1,56 +1,65 @@
-VERSION   1.2
+VERSION	  1.2
 EXPORTS
-	pj_init		@1
-	pj_fwd		@2
-	pj_inv		@3
-	pj_free		@4
-	pj_transform	@5
-        pj_geocentric_to_geodetic @6
-        pj_geodetic_to_geocentric @7
-        pj_deallocate_grids       @8
-        pj_init_plus              @9
-        pj_latlong_from_proj      @10
-        pj_is_latlong             @11
-	pj_get_errno_ref	  @12
-	pj_set_finder             @13
-	pj_strerrno		  @14
-	pj_errno                  @15
-	pj_get_def		  @16
-	pj_dalloc                 @17
-        pj_is_geocent             @18
-        pj_get_release            @19
-	pj_malloc		  @20
-	pj_pr_list		  @21
-	pj_compare_datums	  @22
-	pj_apply_gridshift	  @23
-	pj_datum_transform        @24
-	pj_set_searchpath         @25
-	dmstor			  @26
-	pj_get_ellps_ref          @27
-	pj_get_datums_ref         @28
-	pj_get_units_ref          @29
-	pj_get_list_ref           @30
+	pj_init			@1
+	pj_fwd			@2
+	pj_inv			@3
+	pj_free			@4
+	pj_transform		@5
+	pj_geocentric_to_geodetic @6
+	pj_geodetic_to_geocentric @7
+	pj_deallocate_grids	@8
+	pj_init_plus		@9
+	pj_latlong_from_proj	@10
+	pj_is_latlong		@11
+	pj_get_errno_ref	@12
+	pj_set_finder		@13
+	pj_strerrno		@14
+	pj_errno		@15
+	pj_get_def		@16
+	pj_dalloc		@17
+	pj_is_geocent		@18
+	pj_get_release		@19
+	pj_malloc		@20
+	pj_pr_list		@21
+	pj_compare_datums	@22
+	pj_apply_gridshift	@23
+	pj_datum_transform	@24
+	pj_set_searchpath	@25
+	dmstor			@26
+	pj_get_ellps_ref	@27
+	pj_get_datums_ref	@28
+	pj_get_units_ref	@29
+	pj_get_list_ref		@30
 	pj_get_prime_meridians_ref @31
-	rtodms                    @32
-	set_rtodms                @33
-	pj_factors                @34
-        mk_cheby                  @35
-	adjlon                    @36
-	pj_param		  @37
-	pj_ell_set		  @38
-	pj_mkparam		  @39
-	pj_init_ctx		  @40
-	pj_init_plus_ctx	  @41
-	pj_get_default_ctx        @42
-	pj_get_ctx                @43
-	pj_set_ctx                @44
-	pj_ctx_alloc		  @45
-	pj_ctx_free		  @46
-	pj_ctx_get_errno	  @47
-	pj_ctx_set_errno          @48
-	pj_ctx_set_debug          @49
-	pj_ctx_set_logger         @50
-	pj_ctx_set_app_data       @51
-	pj_ctx_get_app_data       @52
-	pj_log 			  @53
-	pj_clear_initcache @54
+	rtodms			@32
+	set_rtodms		@33
+	pj_factors		@34
+	mk_cheby		@35
+	adjlon			@36
+	pj_param		@37
+	pj_ell_set		@38
+	pj_mkparam		@39
+	pj_init_ctx		@40
+	pj_init_plus_ctx	@41
+	pj_get_default_ctx	@42
+	pj_get_ctx		@43
+	pj_set_ctx		@44
+	pj_ctx_alloc		@45
+	pj_ctx_free		@46
+	pj_ctx_get_errno	@47
+	pj_ctx_set_errno	@48
+	pj_ctx_set_debug	@49
+	pj_ctx_set_logger	@50
+	pj_ctx_set_app_data	@51
+	pj_ctx_get_app_data	@52
+	pj_log			@53
+	pj_clear_initcache	@54
+	geod_init		@55
+	geod_lineinit		@56
+	geod_genposition	@57
+	geod_position		@58
+	geod_gendirect		@59
+	geod_direct		@60
+	geod_geninverse		@61
+	geod_inverse		@62
+	geod_polygonarea	@63