Upgrade geodesic library from GeographicLib 1.46.

* upgrade geodesic.[ch3]
* add test suite geodtest.c and invoke via cmake's add_test
* increment version to 4.9.3 and library version to 11.0.0

1 files changed, 185 insertions, 90 deletions

diff --git a/src/geodesic.c b/src/geodesic.c
index 2dee45e0..8d9c9285 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -14,11 +14,11 @@
  *   Algorithms for geodesics,
  *   J. Geodesy <b>87</b>, 43--55 (2013);
  *   https://dx.doi.org/10.1007/s00190-012-0578-z
- *   Addenda: http://geographiclib.sf.net/geod-addenda.html
+ *   Addenda: http://geographiclib.sourceforge.net/geod-addenda.html
  *
  * See the comments in geodesic.h for documentation.
  *
- * Copyright (c) Charles Karney (2012-2015) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2016) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * http://geographiclib.sourceforge.net/
  */
@@ -119,6 +119,10 @@ static real atanhx(real x) {
   return x < 0 ? -y : y;
 }
 
+static real copysignx(real x, real y) {
+  return fabs(x) * (y < 0 || (y == 0 && 1/y < 0) ? -1 : 1);
+}
+
 static real hypotx(real x, real y)
 { return sqrt(x * x + y * y); }
 
@@ -133,7 +137,7 @@ static real sumx(real u, real v, real* t) {
   volatile real vpp = s - up;
   up -= u;
   vpp -= v;
-  *t = -(up + vpp);
+  if (t) *t = -(up + vpp);
   /* error-free sum:
    * u + v =       s      + t
    *       = round(u + v) + t */
@@ -146,11 +150,12 @@ static real polyval(int N, const real p[], real x) {
   return y;
 }
 
-static real minx(real x, real y)
-{ return x < y ? x : y; }
+/* mimic C++ std::min and std::max */
+static real minx(real a, real b)
+{ return (b < a) ? b : a; }
 
-static real maxx(real x, real y)
-{ return x > y ? x : y; }
+static real maxx(real a, real b)
+{ return (a < b) ? b : a; }
 
 static void swapx(real* x, real* y)
 { real t = *x; *x = *y; *y = t; }
@@ -169,7 +174,7 @@ static real AngNormalize(real x) {
 static real LatFix(real x)
 { return fabs(x) > 90 ? NaN : x; }
 
-static real AngDiff(real x, real y) {
+static real AngDiff(real x, real y, real* e) {
   real t, d = - AngNormalize(sumx(AngNormalize(x), AngNormalize(-y), &t));
   /* Here y - x = d - t (mod 360), exactly, where d is in (-180,180] and
    * abs(t) <= eps (eps = 2^-45 for doubles).  The only case where the
@@ -177,15 +182,16 @@ static real AngDiff(real x, real y) {
    * and t < 0.  The case, d = -180 + eps, t = eps, can't happen, since
    * sum would have returned the exact result in such a case (i.e., given t
    * = 0). */
-  return (d == 180 && t < 0 ? -180 : d) - t;
+  return sumx(d == 180 && t < 0 ? -180 : d, -t, e);
 }
 
 static real AngRound(real x) {
   const real z = 1/(real)(16);
+  if (x == 0) return 0;
   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 ? 0 - y : y;
+  return x < 0 ? -y : y;
 }
 
 static void sincosdx(real x, real* sinx, real* cosx) {
@@ -249,7 +255,7 @@ static real Astroid(real x, real y);
 static real InverseStart(const struct geod_geodesic* g,
                          real sbet1, real cbet1, real dn1,
                          real sbet2, real cbet2, real dn2,
-                         real lam12,
+                         real lam12, real slam12, real clam12,
                          real* psalp1, real* pcalp1,
                          /* Only updated if return val >= 0 */
                          real* psalp2, real* pcalp2,
@@ -261,11 +267,13 @@ static real Lambda12(const struct geod_geodesic* g,
                      real sbet1, real cbet1, real dn1,
                      real sbet2, real cbet2, real dn2,
                      real salp1, real calp1,
+                     real slam120, real clam120,
                      real* psalp2, real* pcalp2,
                      real* psig12,
                      real* pssig1, real* pcsig1,
                      real* pssig2, real* pcsig2,
-                     real* peps, real* pdomg12,
+                     real* peps,
+                     real* psomg12, real* pcomg12,
                      boolx diffp, real* pdlam12,
                      /* Scratch area of the right size */
                      real Ca[]);
@@ -288,7 +296,7 @@ static void accneg(real s[]);
 void geod_init(struct geod_geodesic* g, real a, real f) {
   if (!init) Init();
   g->a = a;
-  g->f = f <= 1 ? f : 1/f;
+  g->f = f;
   g->f1 = 1 - g->f;
   g->e2 = g->f * (2 - g->f);
   g->ep2 = g->e2 / sq(g->f1);   /* e2 / (1 - e2) */
@@ -315,9 +323,11 @@ void geod_init(struct geod_geodesic* g, real a, real f) {
   C4coeff(g);
 }
 
-void geod_lineinit(struct geod_geodesicline* l,
-                   const struct geod_geodesic* g,
-                   real lat1, real lon1, real azi1, unsigned caps) {
+static void geod_lineinit_int(struct geod_geodesicline* l,
+                              const struct geod_geodesic* g,
+                              real lat1, real lon1,
+                              real azi1, real salp1, real calp1,
+                              unsigned caps) {
   real cbet1, sbet1, eps;
   l->a = g->a;
   l->f = g->f;
@@ -331,9 +341,9 @@ void geod_lineinit(struct geod_geodesicline* l,
 
   l->lat1 = LatFix(lat1);
   l->lon1 = lon1;
-  l->azi1 = AngNormalize(azi1);
-  /* Guard against underflow in salp0 */
-  sincosdx(AngRound(l->azi1), &l->salp1, &l->calp1);
+  l->azi1 = azi1;
+  l->salp1 = salp1;
+  l->calp1 = calp1;
 
   sincosdx(AngRound(l->lat1), &sbet1, &cbet1); sbet1 *= l->f1;
   /* Ensure cbet1 = +epsilon at poles */
@@ -396,6 +406,34 @@ void geod_lineinit(struct geod_geodesicline* l,
     l->A4 = sq(l->a) * l->calp0 * l->salp0 * g->e2;
     l->B41 = SinCosSeries(FALSE, l->ssig1, l->csig1, l->C4a, nC4);
   }
+
+  l->a13 = l->s13 = NaN;
+}
+
+void geod_lineinit(struct geod_geodesicline* l,
+                   const struct geod_geodesic* g,
+                   real lat1, real lon1, real azi1, unsigned caps) {
+  azi1 = AngNormalize(azi1);
+  real salp1, calp1;
+  /* Guard against underflow in salp0 */
+  sincosdx(AngRound(azi1), &salp1, &calp1);
+  geod_lineinit_int(l, g, lat1, lon1, azi1, salp1, calp1, caps);
+}
+
+void geod_gendirectline(struct geod_geodesicline* l,
+                        const struct geod_geodesic* g,
+                        real lat1, real lon1, real azi1,
+                        unsigned flags, real a12_s12,
+                        unsigned caps) {
+  geod_lineinit(l, g, lat1, lon1, azi1, caps);
+  geod_gensetdistance(l, flags, a12_s12);
+}
+
+void geod_directline(struct geod_geodesicline* l,
+                        const struct geod_geodesic* g,
+                        real lat1, real lon1, real azi1,
+                        real s12, unsigned caps) {
+  geod_gendirectline(l, g, lat1, lon1, azi1, GEOD_NOFLAGS, s12, caps);
 }
 
 real geod_genposition(const struct geod_geodesicline* l,
@@ -421,7 +459,7 @@ real geod_genposition(const struct geod_geodesicline* l,
 
   outmask &= l->caps & OUT_ALL;
   if (!( TRUE /*Init()*/ &&
-         (flags & GEOD_ARCMODE || (l->caps & GEOD_DISTANCE_IN & OUT_ALL)) ))
+         (flags & GEOD_ARCMODE || (l->caps & (GEOD_DISTANCE_IN & OUT_ALL))) ))
     /* Uninitialized or impossible distance calculation requested */
     return NaN;
 
@@ -464,10 +502,9 @@ real geod_genposition(const struct geod_geodesicline* l,
        *      1/20   5376 146e3  10e3
        *      1/10  829e3  22e6 1.5e6
        *      1/5   157e6 3.8e9 280e6 */
-      real
-        ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12,
-        csig2 = l->csig1 * csig12 - l->ssig1 * ssig12,
-        serr;
+      real serr;
+      ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12;
+      csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;
       B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
       serr = (1 + l->A1m1) * (sig12 + (B12 - l->B11)) - s12_a12 / l->b;
       sig12 = sig12 - serr / sqrt(1 + l->k2 * sq(ssig2));
@@ -499,7 +536,7 @@ real geod_genposition(const struct geod_geodesicline* l,
     s12 = flags & GEOD_ARCMODE ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12;
 
   if (outmask & GEOD_LONGITUDE) {
-    int E = l->salp0 < 0 ? -1 : 1; /* east or west going? */
+    real E = copysignx(1, l->salp0); /* east or west going? */
     /* tan(omg2) = sin(alp0) * tan(sig2) */
     somg2 = l->salp0 * ssig2; comg2 = csig2;  /* No need to normalize */
     /* omg12 = omg2 - omg1 */
@@ -548,14 +585,6 @@ real geod_genposition(const struct geod_geodesicline* l,
       /* alp12 = alp2 - alp1, used in atan2 so no need to normalize */
       salp12 = salp2 * l->calp1 - calp2 * l->salp1;
       calp12 = calp2 * l->calp1 + salp2 * l->salp1;
-      /* The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz
-       * salp12 = -0 and alp12 = -180.  However this depends on the sign being
-       * attached to 0 correctly.  The following ensures the correct
-       * behavior. */
-      if (salp12 == 0 && calp12 < 0) {
-        salp12 = tiny * l->calp1;
-        calp12 = -1;
-      }
     } else {
       /* tan(alp) = tan(alp0) * sec(sig)
        * tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1)
@@ -593,6 +622,21 @@ real geod_genposition(const struct geod_geodesicline* l,
   return flags & GEOD_ARCMODE ? s12_a12 : sig12 / degree;
 }
 
+void geod_setdistance(struct geod_geodesicline* l, real s13) {
+  l->s13 = s13;
+  l->a13 = geod_genposition(l, GEOD_NOFLAGS, l->s13, 0, 0, 0, 0, 0, 0, 0, 0);
+}
+
+static void geod_setarc(struct geod_geodesicline* l, real a13) {
+  l->a13 = a13; l->s13 = NaN;
+  geod_genposition(l, GEOD_ARCMODE, l->a13, 0, 0, 0, &l->s13, 0, 0, 0, 0);
+}
+
+void geod_gensetdistance(struct geod_geodesicline* l,
+ unsigned flags, real s13_a13) {
+  flags & GEOD_ARCMODE ? geod_setarc(l, s13_a13) : geod_setdistance(l, s13_a13);
+}
+
 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);
@@ -630,23 +674,26 @@ void geod_direct(const struct geod_geodesic* g,
                  0, 0, 0, 0, 0);
 }
 
-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;
+static real geod_geninverse_int(const struct geod_geodesic* g,
+                                real lat1, real lon1, real lat2, real lon2,
+                                real* ps12,
+                                real* psalp1, real* pcalp1,
+                                real* psalp2, real* pcalp2,
+                                real* pm12, real* pM12, real* pM21,
+                                real* pS12) {
+  real s12 = 0, m12 = 0, M12 = 0, M21 = 0, S12 = 0;
+  real lon12, lon12s;
   int latsign, lonsign, swapp;
   real sbet1, cbet1, sbet2, cbet2, s12x = 0, m12x = 0;
   real dn1, dn2, lam12, slam12, clam12;
   real a12 = 0, sig12, calp1 = 0, salp1 = 0, calp2 = 0, salp2 = 0;
   real Ca[nC];
   boolx meridian;
-  real omg12 = 0;
+  /* somg12 > 1 marks that it needs to be calculated */
+  real omg12 = 0, somg12 = 2, comg12 = 0;
 
   unsigned outmask =
     (ps12 ? GEOD_DISTANCE : 0U) |
-    (pazi1 || pazi2 ? GEOD_AZIMUTH : 0U) |
     (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
     (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
     (pS12 ? GEOD_AREA : 0U);
@@ -655,16 +702,25 @@ real geod_geninverse(const struct geod_geodesic* g,
   /* 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. */
-  /* If very close to being on the same half-meridian, then make it so. */
-  lon12 = AngRound(AngDiff(lon1, lon2));
+  lon12 = AngDiff(lon1, lon2, &lon12s);
   /* Make longitude difference positive. */
   lonsign = lon12 >= 0 ? 1 : -1;
-  lon12 *= lonsign;
+  /* If very close to being on the same half-meridian, then make it so. */
+  lon12 = lonsign * AngRound(lon12);
+  lon12s = AngRound((180 - lon12) - lonsign * lon12s);
+  lam12 = lon12 * degree;
+  if (lon12 > 90) {
+    sincosdx(lon12s, &slam12, &clam12);
+    clam12 = -clam12;
+  } else
+    sincosdx(lon12, &slam12, &clam12);
+
   /* If really close to the equator, treat as on equator. */
   lat1 = AngRound(LatFix(lat1));
   lat2 = AngRound(LatFix(lat2));
-  /* Swap points so that point with higher (abs) latitude is point 1 */
-  swapp = fabs(lat1) >= fabs(lat2) ? 1 : -1;
+  /* Swap points so that point with higher (abs) latitude is point 1
+   * If one latitude is a nan, then it becomes lat1. */
+  swapp = fabs(lat1) < fabs(lat2) ? -1 : 1;
   if (swapp < 0) {
     lonsign *= -1;
     swapx(&lat1, &lat2);
@@ -712,9 +768,6 @@ real geod_geninverse(const struct geod_geodesic* g,
   dn1 = sqrt(1 + g->ep2 * sq(sbet1));
   dn2 = sqrt(1 + g->ep2 * sq(sbet2));
 
-  lam12 = lon12 * degree;
-  sincosdx(lon12, &slam12, &clam12);
-
   meridian = lat1 == -90 || slam12 == 0;
 
   if (meridian) {
@@ -731,7 +784,7 @@ real geod_geninverse(const struct geod_geodesic* g,
     ssig2 = sbet2; csig2 = calp2 * cbet2;
 
     /* sig12 = sig2 - sig1 */
-    sig12 = atan2(maxx(csig1 * ssig2 - ssig1 * csig2, (real)(0)),
+    sig12 = atan2(maxx((real)(0), csig1 * ssig2 - ssig1 * csig2),
                   csig1 * csig2 + ssig1 * ssig2);
     Lengths(g, g->n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
             cbet1, cbet2, &s12x, &m12x, 0,
@@ -760,7 +813,7 @@ real geod_geninverse(const struct geod_geodesic* g,
   if (!meridian &&
       sbet1 == 0 &&           /* and sbet2 == 0 */
       /* Mimic the way Lambda12 works with calp1 = 0 */
-      (g->f <= 0 || lam12 <= pi - g->f * pi)) {
+      (g->f <= 0 || lon12s >= g->f * 180)) {
 
     /* Geodesic runs along equator */
     calp1 = calp2 = 0; salp1 = salp2 = 1;
@@ -779,7 +832,7 @@ real geod_geninverse(const struct geod_geodesic* g,
     /* Figure a starting point for Newton's method */
     real dnm = 0;
     sig12 = InverseStart(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
-                         lam12,
+                         lam12, slam12, clam12,
                          &salp1, &calp1, &salp2, &calp2, &dnm,
                          Ca);
 
@@ -813,13 +866,13 @@ real geod_geninverse(const struct geod_geodesic* g,
         /* the WGS84 test set: mean = 1.47, sd = 1.25, max = 16
          * WGS84 and random input: mean = 2.85, sd = 0.60 */
         real dv = 0,
-          v = (Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
+          v = Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
+                        slam12, clam12,
                         &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2,
-                        &eps, &omg12, numit < maxit1, &dv, Ca)
-               - lam12);
+                        &eps, &somg12, &comg12, numit < maxit1, &dv, Ca);
         /* 2 * tol0 is approximately 1 ulp for a number in [0, pi]. */
         /* Reversed test to allow escape with NaNs */
-        if (tripb || !(fabs(v) >= (tripn ? 8 : 2) * tol0)) break;
+        if (tripb || !(fabs(v) >= (tripn ? 8 : 1) * tol0)) break;
         /* Update bracketing values */
         if (v > 0 && (numit > maxit1 || calp1/salp1 > calp1b/salp1b))
           { salp1b = salp1; calp1b = calp1; }
@@ -864,7 +917,6 @@ real geod_geninverse(const struct geod_geodesic* g,
       m12x *= g->b;
       s12x *= g->b;
       a12 = sig12 / degree;
-      omg12 = lam12 - omg12;
     }
   }
 
@@ -900,15 +952,22 @@ real geod_geninverse(const struct geod_geodesic* g,
       /* Avoid problems with indeterminate sig1, sig2 on equator */
       S12 = 0;
 
+    if (!meridian) {
+      if (somg12 > 1) {
+        somg12 = sin(omg12); comg12 = cos(omg12);
+      } else
+        norm2(&somg12, &comg12);
+    }
+
     if (!meridian &&
-        omg12 < (real)(0.75) * pi &&   /* Long difference too big */
-        sbet2 - sbet1 < (real)(1.75)) { /* Lat difference too big */
+        /* omg12 < 3/4 * pi */
+        comg12 > -(real)(0.7071) &&     /* Long difference not too big */
+        sbet2 - sbet1 < (real)(1.75)) { /* Lat difference not too big */
       /* Use tan(Gamma/2) = tan(omg12/2)
        * * (tan(bet1/2)+tan(bet2/2))/(1+tan(bet1/2)*tan(bet2/2))
        * with tan(x/2) = sin(x)/(1+cos(x)) */
       real
-        somg12 = sin(omg12), domg12 = 1 + cos(omg12),
-        dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
+        domg12 = 1 + comg12, dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
       alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
                          domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
     } else {
@@ -943,18 +1002,13 @@ real geod_geninverse(const struct geod_geodesic* g,
   salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
   salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
 
-  if (outmask & GEOD_AZIMUTH) {
-    /* minus signs give range [-180, 180). 0- converts -0 to +0. */
-    azi1 = atan2dx(salp1, calp1);
-    azi2 = atan2dx(salp2, calp2);
-  }
+  if (psalp1) *psalp1 = salp1;
+  if (pcalp1) *pcalp1 = calp1;
+  if (psalp2) *psalp2 = salp2;
+  if (pcalp2) *pcalp2 = calp2;
 
   if (outmask & GEOD_DISTANCE)
     *ps12 = s12;
-  if (outmask & GEOD_AZIMUTH) {
-    if (pazi1) *pazi1 = azi1;
-    if (pazi2) *pazi2 = azi2;
-  }
   if (outmask & GEOD_REDUCEDLENGTH)
     *pm12 = m12;
   if (outmask & GEOD_GEODESICSCALE) {
@@ -968,6 +1022,35 @@ real geod_geninverse(const struct geod_geodesic* g,
   return a12;
 }
 
+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 salp1, calp1, salp2, calp2,
+    a12 = geod_geninverse_int(g, lat1, lon1, lat2, lon2, ps12,
+                              &salp1, &calp1, &salp2, &calp2,
+                              pm12, pM12, pM21, pS12);
+  if (pazi1) *pazi1 = atan2dx(salp1, calp1);
+  if (pazi2) *pazi2 = atan2dx(salp2, calp2);
+  return a12;
+}
+
+void geod_inverseline(struct geod_geodesicline* l,
+                      const struct geod_geodesic* g,
+                      real lat1, real lon1, real lat2, real lon2,
+                      unsigned caps) {
+  real salp1, calp1,
+    a12 = geod_geninverse_int(g, lat1, lon1, lat2, lon2, 0,
+                              &salp1, &calp1, 0, 0,
+                              0, 0, 0, 0),
+    azi1 = atan2dx(salp1, calp1);
+  caps = caps ? caps : GEOD_DISTANCE_IN | GEOD_LONGITUDE;
+  /* Ensure that a12 can be converted to a distance */
+  if (caps & (OUT_ALL & GEOD_DISTANCE_IN)) caps |= GEOD_DISTANCE;
+  geod_lineinit_int(l, g, lat1, lon1, azi1, salp1, calp1, caps);
+  geod_setarc(l, a12);
+}
+
 void geod_inverse(const struct geod_geodesic* g,
                   real lat1, real lon1, real lat2, real lon2,
                   real* ps12, real* pazi1, real* pazi2) {
@@ -1112,7 +1195,7 @@ real Astroid(real x, real y) {
 real InverseStart(const struct geod_geodesic* g,
                   real sbet1, real cbet1, real dn1,
                   real sbet2, real cbet2, real dn2,
-                  real lam12,
+                  real lam12, real slam12, real clam12,
                   real* psalp1, real* pcalp1,
                   /* Only updated if return val >= 0 */
                   real* psalp2, real* pcalp2,
@@ -1133,7 +1216,7 @@ real InverseStart(const struct geod_geodesic* g,
   real sbet12a;
   boolx shortline = cbet12 >= 0 && sbet12 < (real)(0.5) &&
     cbet2 * lam12 < (real)(0.5);
-  real omg12 = lam12, somg12, comg12, ssig12, csig12;
+  real somg12, comg12, ssig12, csig12;
 #if defined(__GNUC__) && __GNUC__ == 4 &&       \
   (__GNUC_MINOR__ < 6 || defined(__MINGW32__))
   /* Volatile declaration needed to fix inverse cases
@@ -1151,14 +1234,16 @@ real InverseStart(const struct geod_geodesic* g,
   sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
 #endif
   if (shortline) {
-    real sbetm2 = sq(sbet1 + sbet2);
+    real sbetm2 = sq(sbet1 + sbet2), omg12;
     /* sin((bet1+bet2)/2)^2
      * =  (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2) */
     sbetm2 /= sbetm2 + sq(cbet1 + cbet2);
     dnm = sqrt(1 + g->ep2 * sbetm2);
-    omg12 /= g->f1 * dnm;
+    omg12 = lam12 / (g->f1 * dnm);
+    somg12 = sin(omg12); comg12 = cos(omg12);
+  } else {
+    somg12 = slam12; comg12 = clam12;
   }
-  somg12 = sin(omg12); comg12 = cos(omg12);
 
   salp1 = cbet2 * somg12;
   calp1 = comg12 >= 0 ?
@@ -1188,6 +1273,7 @@ real InverseStart(const struct geod_geodesic* g,
      * 56.320923501171 0 -56.320923501171 179.664747671772880215
      * which otherwise fails with g++ 4.4.4 x86 -O3 */
     volatile real x;
+    real lam12x = atan2(-slam12, -clam12); /* lam12 - pi */
     if (g->f >= 0) {            /* In fact f == 0 does not get here */
       /* x = dlong, y = dlat */
       {
@@ -1198,7 +1284,7 @@ real InverseStart(const struct geod_geodesic* g,
       }
       betscale = lamscale * cbet1;
 
-      x = (lam12 - pi) / lamscale;
+      x = lam12x / lamscale;
       y = sbet12a / betscale;
     } else {                    /* f < 0 */
       /* x = dlat, y = dlong */
@@ -1215,7 +1301,7 @@ real InverseStart(const struct geod_geodesic* g,
       betscale = x < -(real)(0.01) ? sbet12a / x :
         -g->f * sq(cbet1) * pi;
       lamscale = betscale / cbet1;
-      y = (lam12 - pi) / lamscale;
+      y = lam12x / lamscale;
     }
 
     if (y > -tol1 && x > -1 - xthresh) {
@@ -1292,19 +1378,22 @@ real Lambda12(const struct geod_geodesic* g,
               real sbet1, real cbet1, real dn1,
               real sbet2, real cbet2, real dn2,
               real salp1, real calp1,
+              real slam120, real clam120,
               real* psalp2, real* pcalp2,
               real* psig12,
               real* pssig1, real* pcsig1,
               real* pssig2, real* pcsig2,
-              real* peps, real* pdomg12,
+              real* peps,
+              real* psomg12, real* pcomg12,
               boolx diffp, real* pdlam12,
               /* Scratch area of the right size */
               real Ca[]) {
   real salp2 = 0, calp2 = 0, sig12 = 0,
-    ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0, domg12 = 0, dlam12 = 0;
+    ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0,
+    somg12 = 0, comg12 = 0, dlam12 = 0;
   real salp0, calp0;
-  real somg1, comg1, somg2, comg2, omg12, lam12;
-  real B312, h0, k2;
+  real somg1, comg1, somg2, comg2, lam12;
+  real B312, eta, k2;
 
   if (sbet1 == 0 && calp1 == 0)
     /* Break degeneracy of equatorial line.  This case has already been
@@ -1345,20 +1434,21 @@ real Lambda12(const struct geod_geodesic* g,
   /* norm2(&somg2, &comg2); -- don't need to normalize! */
 
   /* sig12 = sig2 - sig1, limit to [0, pi] */
-  sig12 = atan2(maxx(csig1 * ssig2 - ssig1 * csig2, (real)(0)),
+  sig12 = atan2(maxx((real)(0), csig1 * ssig2 - ssig1 * csig2),
                 csig1 * csig2 + ssig1 * ssig2);
 
   /* omg12 = omg2 - omg1, limit to [0, pi] */
-  omg12 = atan2(maxx(comg1 * somg2 - somg1 * comg2, (real)(0)),
-                comg1 * comg2 + somg1 * somg2);
+  somg12 = maxx((real)(0), comg1 * somg2 - somg1 * comg2);
+  comg12 =                 comg1 * comg2 + somg1 * somg2;
+  /* eta = omg12 - lam120 */
+  eta = atan2(somg12 * clam120 - comg12 * slam120,
+              comg12 * clam120 + somg12 * slam120);
   k2 = sq(calp0) * g->ep2;
   eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
   C3f(g, eps, Ca);
   B312 = (SinCosSeries(TRUE, ssig2, csig2, Ca, nC3-1) -
           SinCosSeries(TRUE, ssig1, csig1, Ca, nC3-1));
-  h0 = -g->f * A3f(g, eps);
-  domg12 = salp0 * h0 * (sig12 + B312);
-  lam12 = omg12 + domg12;
+  lam12 = eta - g->f * A3f(g, eps) * salp0 * (sig12 + B312);
 
   if (diffp) {
     if (calp2 == 0)
@@ -1378,7 +1468,8 @@ real Lambda12(const struct geod_geodesic* g,
   *pssig2 = ssig2;
   *pcsig2 = csig2;
   *peps = eps;
-  *pdomg12 = domg12;
+  *psomg12 = somg12;
+  *pcomg12 = comg12;
   if (diffp)
     *pdlam12 = dlam12;
 
@@ -1653,7 +1744,7 @@ int transit(real lon1, real lon2) {
   /* Compute lon12 the same way as Geodesic::Inverse. */
   lon1 = AngNormalize(lon1);
   lon2 = AngNormalize(lon2);
-  lon12 = AngDiff(lon1, lon2);
+  lon12 = AngDiff(lon1, lon2, 0);
   return lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 :
     (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0);
 }
@@ -1699,8 +1790,12 @@ void accneg(real s[]) {
 }
 
 void geod_polygon_init(struct geod_polygon* p, boolx polylinep) {
-  p->lat0 = p->lon0 = p->lat = p->lon = NaN;
   p->polyline = (polylinep != 0);
+  geod_polygon_clear(p);
+}
+
+void geod_polygon_clear(struct geod_polygon* p) {
+  p->lat0 = p->lon0 = p->lat = p->lon = NaN;
   accini(p->P);
   accini(p->A);
   p->num = p->crossings = 0;