Update to version 1.43 of the geodesic routines. This fixes two

relatively obscure problems.  (1) The business of returning an unrolled
longitude with the solution to the direct problem was broken for
west-going geodesics.  (2) For flattening > 1/100, a slight inaccurate
result was returned for a12 in the direct calculation.


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

1 files changed, 239 insertions, 141 deletions

diff --git a/src/geodesic.c b/src/geodesic.c
index fd0214c7..9a8f043c 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -18,7 +18,7 @@
  *
  * See the comments in geodesic.h for documentation.
  *
- * Copyright (c) Charles Karney (2012-2014) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2015) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * http://geographiclib.sourceforge.net/
  */
@@ -27,8 +27,10 @@
 #include <math.h>
 
 #define GEOGRAPHICLIB_GEODESIC_ORDER 6
+#define nA1   GEOGRAPHICLIB_GEODESIC_ORDER
 #define nC1   GEOGRAPHICLIB_GEODESIC_ORDER
 #define nC1p  GEOGRAPHICLIB_GEODESIC_ORDER
+#define nA2   GEOGRAPHICLIB_GEODESIC_ORDER
 #define nC2   GEOGRAPHICLIB_GEODESIC_ORDER
 #define nA3   GEOGRAPHICLIB_GEODESIC_ORDER
 #define nA3x  nA3
@@ -137,6 +139,12 @@ static real sumx(real u, real v, real* t) {
   return s;
 }
 
+static real polyval(int N, const real p[], real x) {
+  real y = N < 0 ? 0 : *p++;
+  while (--N >= 0) y = y * x + *p++;
+  return y;
+}
+
 static real minx(real x, real y)
 { return x < y ? x : y; }
 
@@ -146,7 +154,7 @@ static real maxx(real x, real y)
 static void swapx(real* x, real* y)
 { real t = *x; *x = *y; *y = t; }
 
-static void SinCosNorm(real* sinx, real* cosx) {
+static void norm2(real* sinx, real* cosx) {
   real r = hypotx(*sinx, *cosx);
   *sinx /= r;
   *cosx /= r;
@@ -171,7 +179,7 @@ static real AngRound(real x) {
   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;
+  return x < 0 ? 0 - y : y;
 }
 
 static void A3coeff(struct geod_geodesic* g);
@@ -270,7 +278,8 @@ void geod_lineinit(struct geod_geodesicline* l,
   l->f1 = g->f1;
   /* If caps is 0 assume the standard direct calculation */
   l->caps = (caps ? caps : GEOD_DISTANCE_IN | GEOD_LONGITUDE) |
-    GEOD_LATITUDE | GEOD_AZIMUTH; /* Always allow latitude and azimuth */
+    /* always allow latitude and azimuth and unrolling of longitude */
+    GEOD_LATITUDE | GEOD_AZIMUTH | GEOD_LONG_UNROLL;
 
   l->lat1 = lat1;
   l->lon1 = lon1;
@@ -286,7 +295,7 @@ void geod_lineinit(struct geod_geodesicline* l,
   /* Ensure cbet1 = +epsilon at poles */
   sbet1 = l->f1 * sin(phi);
   cbet1 = fabs(lat1) == 90 ? tiny : cos(phi);
-  SinCosNorm(&sbet1, &cbet1);
+  norm2(&sbet1, &cbet1);
   l->dn1 = sqrt(1 + g->ep2 * sq(sbet1));
 
   /* Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0), */
@@ -305,8 +314,8 @@ void geod_lineinit(struct geod_geodesicline* l,
    * With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi. */
   l->ssig1 = sbet1; l->somg1 = l->salp0 * sbet1;
   l->csig1 = l->comg1 = sbet1 != 0 || l->calp1 != 0 ? cbet1 * l->calp1 : 1;
-  SinCosNorm(&l->ssig1, &l->csig1); /* sig1 in (-pi, pi] */
-  /* SinCosNorm(somg1, comg1); -- don't need to normalize! */
+  norm2(&l->ssig1, &l->csig1); /* sig1 in (-pi, pi] */
+  /* norm2(somg1, comg1); -- don't need to normalize! */
 
   l->k2 = sq(l->calp0) * g->ep2;
   eps = l->k2 / (2 * (1 + sqrt(1 + l->k2)) + l->k2);
@@ -452,12 +461,14 @@ 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? */
     /* tan(omg2) = sin(alp0) * tan(sig2) */
     somg2 = l->salp0 * ssig2; comg2 = csig2;  /* No need to normalize */
     /* omg12 = omg2 - omg1 */
-    omg12 = flags & GEOD_LONG_NOWRAP ? sig12
-      - (atan2(ssig2, csig2) - atan2(l->ssig1, l->csig1))
-      + (atan2(somg2, comg2) - atan2(l->somg1, l->comg1))
+    omg12 = flags & GEOD_LONG_UNROLL
+      ? E * (sig12
+             - (atan2(    ssig2, csig2) - atan2(    l->ssig1, l->csig1))
+             + (atan2(E * somg2, comg2) - atan2(E * l->somg1, l->comg1)))
       : atan2(somg2 * l->comg1 - comg2 * l->somg1,
               comg2 * l->comg1 + somg2 * l->somg1);
     lam12 = omg12 + l->A3c *
@@ -466,7 +477,7 @@ real geod_genposition(const struct geod_geodesicline* l,
     lon12 = lam12 / degree;
     /* Use AngNormalize2 because longitude might have wrapped multiple
      * times. */
-    lon2 = flags & GEOD_LONG_NOWRAP ? l->lon1 + lon12 :
+    lon2 = flags & GEOD_LONG_UNROLL ? l->lon1 + lon12 :
       AngNormalize(AngNormalize(l->lon1) + AngNormalize2(lon12));
   }
 
@@ -499,7 +510,7 @@ real geod_genposition(const struct geod_geodesicline* l,
       B42 = SinCosSeries(FALSE, ssig2, csig2, l->C4a, nC4);
     real salp12, calp12;
     if (l->calp0 == 0 || l->salp0 == 0) {
-      /* alp12 = alp2 - alp1, used in atan2 so no need to normalized */
+      /* 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
@@ -645,13 +656,13 @@ real geod_geninverse(const struct geod_geodesic* g,
   /* Ensure cbet1 = +epsilon at poles */
   sbet1 = g->f1 * sin(phi);
   cbet1 = lat1 == -90 ? tiny : cos(phi);
-  SinCosNorm(&sbet1, &cbet1);
+  norm2(&sbet1, &cbet1);
 
   phi = lat2 * degree;
   /* Ensure cbet2 = +epsilon at poles */
   sbet2 = g->f1 * sin(phi);
   cbet2 = fabs(lat2) == 90 ? tiny : cos(phi);
-  SinCosNorm(&sbet2, &cbet2);
+  norm2(&sbet2, &cbet2);
 
   /* If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the
    * |bet1| - |bet2|.  Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is
@@ -771,7 +782,7 @@ real geod_geninverse(const struct geod_geodesic* g,
       for (tripn = FALSE, tripb = FALSE; numit < maxit2; ++numit) {
         /* the WGS84 test set: mean = 1.47, sd = 1.25, max = 16
          * WGS84 and random input: mean = 2.85, sd = 0.60 */
-        real dv,
+        real dv = 0,
           v = (Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
                         &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2,
                         &eps, &omg12, numit < maxit1, &dv, C1a, C2a, C3a)
@@ -793,7 +804,7 @@ real geod_geninverse(const struct geod_geodesic* g,
           if (nsalp1 > 0 && fabs(dalp1) < pi) {
             calp1 = calp1 * cdalp1 - salp1 * sdalp1;
             salp1 = nsalp1;
-            SinCosNorm(&salp1, &calp1);
+            norm2(&salp1, &calp1);
             /* In some regimes we don't get quadratic convergence because
              * slope -> 0.  So use convergence conditions based on epsilon
              * instead of sqrt(epsilon). */
@@ -811,7 +822,7 @@ real geod_geninverse(const struct geod_geodesic* g,
          * WGS84 and random input: mean = 4.74, sd = 0.99 */
         salp1 = (salp1a + salp1b)/2;
         calp1 = (calp1a + calp1b)/2;
-        SinCosNorm(&salp1, &calp1);
+        norm2(&salp1, &calp1);
         tripn = FALSE;
         tripb = (fabs(salp1a - salp1) + (calp1a - calp1) < tolb ||
                  fabs(salp1 - salp1b) + (calp1 - calp1b) < tolb);
@@ -852,8 +863,8 @@ real geod_geninverse(const struct geod_geodesic* g,
         A4 = sq(g->a) * calp0 * salp0 * g->e2;
       real C4a[nC4];
       real B41, B42;
-      SinCosNorm(&ssig1, &csig1);
-      SinCosNorm(&ssig2, &csig2);
+      norm2(&ssig1, &csig1);
+      norm2(&ssig2, &csig2);
       C4f(g, eps, C4a);
       B41 = SinCosSeries(FALSE, ssig1, csig1, C4a, nC4);
       B42 = SinCosSeries(FALSE, ssig2, csig2, C4a, nC4);
@@ -1119,7 +1130,7 @@ real InverseStart(const struct geod_geodesic* g,
     salp2 = cbet1 * somg12;
     calp2 = sbet12 - cbet1 * sbet2 *
       (comg12 >= 0 ? sq(somg12) / (1 + comg12) : 1 - comg12);
-    SinCosNorm(&salp2, &calp2);
+    norm2(&salp2, &calp2);
     /* Set return value */
     sig12 = atan2(ssig12, csig12);
   } else if (fabs(g->n) > (real)(0.1) || /* No astroid calc if too eccentric */
@@ -1219,7 +1230,7 @@ real InverseStart(const struct geod_geodesic* g,
   }
   /* Sanity check on starting guess.  Backwards check allows NaN through. */
   if (!(salp1 <= 0))
-    SinCosNorm(&salp1, &calp1);
+    norm2(&salp1, &calp1);
   else {
     salp1 = 1; calp1 = 0;
   }
@@ -1266,8 +1277,8 @@ real Lambda12(const struct geod_geodesic* g,
    * tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1) */
   ssig1 = sbet1; somg1 = salp0 * sbet1;
   csig1 = comg1 = calp1 * cbet1;
-  SinCosNorm(&ssig1, &csig1);
-  /* SinCosNorm(&somg1, &comg1); -- don't need to normalize! */
+  norm2(&ssig1, &csig1);
+  /* norm2(&somg1, &comg1); -- don't need to normalize! */
 
   /* Enforce symmetries in the case abs(bet2) = -bet1.  Need to be careful
    * about this case, since this can yield singularities in the Newton
@@ -1288,8 +1299,8 @@ real Lambda12(const struct geod_geodesic* g,
    * tan(omg2) = sin(alp0) * tan(sig2). */
   ssig2 = sbet2; somg2 = salp0 * sbet2;
   csig2 = comg2 = calp2 * cbet2;
-  SinCosNorm(&ssig2, &csig2);
-  /* SinCosNorm(&somg2, &comg2); -- don't need to normalize! */
+  norm2(&ssig2, &csig2);
+  /* norm2(&somg2, &comg2); -- don't need to normalize! */
 
   /* sig12 = sig2 - sig1, limit to [0, pi] */
   sig12 = atan2(maxx(csig1 * ssig2 - ssig1 * csig2, (real)(0)),
@@ -1335,177 +1346,264 @@ real Lambda12(const struct geod_geodesic* g,
 }
 
 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;
-  for (i = nA3x; i; )
-    v = eps * v + g->A3x[--i];
-  return v;
+  /* Evaluate A3 */
+  return polyval(nA3 - 1, g->A3x, eps);
 }
 
 void C3f(const struct geod_geodesic* g, real eps, real c[]) {
-  /* Evaluate C3 coeffs by Horner's method
+  /* Evaluate C3 coeffs
    * Elements c[1] thru c[nC3 - 1] are set */
-  int i, j, k;
   real mult = 1;
-  for (j = nC3x, k = nC3 - 1; k; ) {
-    real t = 0;
-    for (i = nC3 - k; i; --i)
-      t = eps * t + g->C3x[--j];
-    c[k--] = t;
-  }
-
-  for (k = 1; k < nC3; ) {
+  int o = 0, l;
+  for (l = 1; l < nC3; ++l) {   /* l is index of C3[l] */
+    int m = nC3 - l - 1;        /* order of polynomial in eps */
     mult *= eps;
-    c[k++] *= mult;
+    c[l] = mult * polyval(m, g->C3x + o, eps);
+    o += m + 1;
   }
 }
 
 void C4f(const struct geod_geodesic* g, real eps, real c[]) {
-  /* Evaluate C4 coeffs by Horner's method
+  /* Evaluate C4 coeffs
    * Elements c[0] thru c[nC4 - 1] are set */
-  int i, j, k;
   real mult = 1;
-  for (j = nC4x, k = nC4; k; ) {
-    real t = 0;
-    for (i = nC4 - k + 1; i; --i)
-      t = eps * t + g->C4x[--j];
-    c[--k] = t;
-  }
-
-  for (k = 1; k < nC4; ) {
+  int o = 0, l;
+  for (l = 0; l < nC4; ++l) {   /* l is index of C4[l] */
+    int m = nC4 - l - 1;        /* order of polynomial in eps */
+    c[l] = mult * polyval(m, g->C4x + o, eps);
+    o += m + 1;
     mult *= eps;
-    c[k++] *= mult;
   }
 }
 
-/* Generated by Maxima on 2010-09-04 10:26:17-04:00 */
-
 /* The scale factor A1-1 = mean value of (d/dsigma)I1 - 1 */
 real A1m1f(real eps)  {
-  real
-    eps2 = sq(eps),
-    t = eps2*(eps2*(eps2+4)+64)/256;
+  static const real coeff[] = {
+    /* (1-eps)*A1-1, polynomial in eps2 of order 3 */
+    1, 4, 64, 0, 256,
+  };
+  int m = nA1/2;
+  real t = polyval(m, coeff, sq(eps)) / coeff[m + 1];
   return (t + eps) / (1 - eps);
 }
 
 /* The coefficients C1[l] in the Fourier expansion of B1 */
 void C1f(real eps, real c[])  {
+  static const real coeff[] = {
+    /* C1[1]/eps^1, polynomial in eps2 of order 2 */
+    -1, 6, -16, 32,
+    /* C1[2]/eps^2, polynomial in eps2 of order 2 */
+    -9, 64, -128, 2048,
+    /* C1[3]/eps^3, polynomial in eps2 of order 1 */
+    9, -16, 768,
+    /* C1[4]/eps^4, polynomial in eps2 of order 1 */
+    3, -5, 512,
+    /* C1[5]/eps^5, polynomial in eps2 of order 0 */
+    -7, 1280,
+    /* C1[6]/eps^6, polynomial in eps2 of order 0 */
+    -7, 2048,
+  };
   real
     eps2 = sq(eps),
     d = eps;
-  c[1] = d*((6-eps2)*eps2-16)/32;
-  d *= eps;
-  c[2] = d*((64-9*eps2)*eps2-128)/2048;
-  d *= eps;
-  c[3] = d*(9*eps2-16)/768;
-  d *= eps;
-  c[4] = d*(3*eps2-5)/512;
-  d *= eps;
-  c[5] = -7*d/1280;
-  d *= eps;
-  c[6] = -7*d/2048;
+  int o = 0, l;
+  for (l = 1; l <= nC1; ++l) {  /* l is index of C1p[l] */
+    int m = (nC1 - l) / 2;      /* order of polynomial in eps^2 */
+    c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
+    o += m + 2;
+    d *= eps;
+  }
 }
 
 /* The coefficients C1p[l] in the Fourier expansion of B1p */
 void C1pf(real eps, real c[])  {
+  static const real coeff[] = {
+    /* C1p[1]/eps^1, polynomial in eps2 of order 2 */
+    205, -432, 768, 1536,
+    /* C1p[2]/eps^2, polynomial in eps2 of order 2 */
+    4005, -4736, 3840, 12288,
+    /* C1p[3]/eps^3, polynomial in eps2 of order 1 */
+    -225, 116, 384,
+    /* C1p[4]/eps^4, polynomial in eps2 of order 1 */
+    -7173, 2695, 7680,
+    /* C1p[5]/eps^5, polynomial in eps2 of order 0 */
+    3467, 7680,
+    /* C1p[6]/eps^6, polynomial in eps2 of order 0 */
+    38081, 61440,
+  };
   real
     eps2 = sq(eps),
     d = eps;
-  c[1] = d*(eps2*(205*eps2-432)+768)/1536;
-  d *= eps;
-  c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288;
-  d *= eps;
-  c[3] = d*(116-225*eps2)/384;
-  d *= eps;
-  c[4] = d*(2695-7173*eps2)/7680;
-  d *= eps;
-  c[5] = 3467*d/7680;
-  d *= eps;
-  c[6] = 38081*d/61440;
+  int o = 0, l;
+  for (l = 1; l <= nC1p; ++l) { /* l is index of C1p[l] */
+    int m = (nC1p - l) / 2;     /* order of polynomial in eps^2 */
+    c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
+    o += m + 2;
+    d *= eps;
+  }
 }
 
 /* The scale factor A2-1 = mean value of (d/dsigma)I2 - 1 */
 real A2m1f(real eps)  {
-  real
-    eps2 = sq(eps),
-    t = eps2*(eps2*(25*eps2+36)+64)/256;
+  static const real coeff[] = {
+    /* A2/(1-eps)-1, polynomial in eps2 of order 3 */
+    25, 36, 64, 0, 256,
+  };
+  int m = nA2/2;
+  real t = polyval(m, coeff, sq(eps)) / coeff[m + 1];
   return t * (1 - eps) - eps;
 }
 
 /* The coefficients C2[l] in the Fourier expansion of B2 */
 void C2f(real eps, real c[])  {
+  static const real coeff[] = {
+    /* C2[1]/eps^1, polynomial in eps2 of order 2 */
+    1, 2, 16, 32,
+    /* C2[2]/eps^2, polynomial in eps2 of order 2 */
+    35, 64, 384, 2048,
+    /* C2[3]/eps^3, polynomial in eps2 of order 1 */
+    15, 80, 768,
+    /* C2[4]/eps^4, polynomial in eps2 of order 1 */
+    7, 35, 512,
+    /* C2[5]/eps^5, polynomial in eps2 of order 0 */
+    63, 1280,
+    /* C2[6]/eps^6, polynomial in eps2 of order 0 */
+    77, 2048,
+  };
   real
     eps2 = sq(eps),
     d = eps;
-  c[1] = d*(eps2*(eps2+2)+16)/32;
-  d *= eps;
-  c[2] = d*(eps2*(35*eps2+64)+384)/2048;
-  d *= eps;
-  c[3] = d*(15*eps2+80)/768;
-  d *= eps;
-  c[4] = d*(7*eps2+35)/512;
-  d *= eps;
-  c[5] = 63*d/1280;
-  d *= eps;
-  c[6] = 77*d/2048;
+  int o = 0, l;
+  for (l = 1; l <= nC2; ++l) { /* l is index of C2[l] */
+    int m = (nC2 - l) / 2;     /* order of polynomial in eps^2 */
+    c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
+    o += m + 2;
+    d *= eps;
+  }
 }
 
 /* The scale factor A3 = mean value of (d/dsigma)I3 */
 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;
-  g->A3x[3] = ((-g->n-3)*g->n-1)/16;
-  g->A3x[4] = (-2*g->n-3)/64;
-  g->A3x[5] = -3/(real)(128);
+  static const real coeff[] = {
+    /* A3, coeff of eps^5, polynomial in n of order 0 */
+    -3, 128,
+    /* A3, coeff of eps^4, polynomial in n of order 1 */
+    -2, -3, 64,
+    /* A3, coeff of eps^3, polynomial in n of order 2 */
+    -1, -3, -1, 16,
+    /* A3, coeff of eps^2, polynomial in n of order 2 */
+    3, -1, -2, 8,
+    /* A3, coeff of eps^1, polynomial in n of order 1 */
+    1, -1, 2,
+    /* A3, coeff of eps^0, polynomial in n of order 0 */
+    1, 1,
+  };
+  int o = 0, k = 0, j;
+  for (j = nA3 - 1; j >= 0; --j) {             /* coeff of eps^j */
+    int m = nA3 - j - 1 < j ? nA3 - j - 1 : j; /* order of polynomial in n */
+    g->A3x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
+    o += m + 2;
+  }
 }
 
 /* The coefficients C3[l] in the Fourier expansion of B3 */
 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;
-  g->C3x[3] = (2*g->n+5)/128;
-  g->C3x[4] = 3/(real)(128);
-  g->C3x[5] = ((g->n-3)*g->n+2)/32;
-  g->C3x[6] = ((-3*g->n-2)*g->n+3)/64;
-  g->C3x[7] = (g->n+3)/128;
-  g->C3x[8] = 5/(real)(256);
-  g->C3x[9] = (g->n*(5*g->n-9)+5)/192;
-  g->C3x[10] = (9-10*g->n)/384;
-  g->C3x[11] = 7/(real)(512);
-  g->C3x[12] = (7-14*g->n)/512;
-  g->C3x[13] = 7/(real)(512);
-  g->C3x[14] = 21/(real)(2560);
+  static const real coeff[] = {
+    /* C3[1], coeff of eps^5, polynomial in n of order 0 */
+    3, 128,
+    /* C3[1], coeff of eps^4, polynomial in n of order 1 */
+    2, 5, 128,
+    /* C3[1], coeff of eps^3, polynomial in n of order 2 */
+    -1, 3, 3, 64,
+    /* C3[1], coeff of eps^2, polynomial in n of order 2 */
+    -1, 0, 1, 8,
+    /* C3[1], coeff of eps^1, polynomial in n of order 1 */
+    -1, 1, 4,
+    /* C3[2], coeff of eps^5, polynomial in n of order 0 */
+    5, 256,
+    /* C3[2], coeff of eps^4, polynomial in n of order 1 */
+    1, 3, 128,
+    /* C3[2], coeff of eps^3, polynomial in n of order 2 */
+    -3, -2, 3, 64,
+    /* C3[2], coeff of eps^2, polynomial in n of order 2 */
+    1, -3, 2, 32,
+    /* C3[3], coeff of eps^5, polynomial in n of order 0 */
+    7, 512,
+    /* C3[3], coeff of eps^4, polynomial in n of order 1 */
+    -10, 9, 384,
+    /* C3[3], coeff of eps^3, polynomial in n of order 2 */
+    5, -9, 5, 192,
+    /* C3[4], coeff of eps^5, polynomial in n of order 0 */
+    7, 512,
+    /* C3[4], coeff of eps^4, polynomial in n of order 1 */
+    -14, 7, 512,
+    /* C3[5], coeff of eps^5, polynomial in n of order 0 */
+    21, 2560,
+  };
+  int o = 0, k = 0, l, j;
+  for (l = 1; l < nC3; ++l) {                    /* l is index of C3[l] */
+    for (j = nC3 - 1; j >= l; --j) {             /* coeff of eps^j */
+      int m = nC3 - j - 1 < j ? nC3 - j - 1 : j; /* order of polynomial in n */
+      g->C3x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
+      o += m + 2;
+    }
+  }
 }
 
-/* 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 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;
-  g->C4x[2] = (g->n*((14144-10656*g->n)*g->n-4576)-858)/45045;
-  g->C4x[3] = ((-224*g->n-4784)*g->n+1573)/45045;
-  g->C4x[4] = (1088*g->n+156)/45045;
-  g->C4x[5] = 97/(real)(15015);
-  g->C4x[6] = (g->n*(g->n*((-64*g->n-624)*g->n+4576)-6864)+3003)/135135;
-  g->C4x[7] = (g->n*(g->n*(5952*g->n-11648)+9152)-2574)/135135;
-  g->C4x[8] = (g->n*(5792*g->n+1040)-1287)/135135;
-  g->C4x[9] = (468-2944*g->n)/135135;
-  g->C4x[10] = 1/(real)(9009);
-  g->C4x[11] = (g->n*((4160-1440*g->n)*g->n-4576)+1716)/225225;
-  g->C4x[12] = ((4992-8448*g->n)*g->n-1144)/225225;
-  g->C4x[13] = (1856*g->n-936)/225225;
-  g->C4x[14] = 8/(real)(10725);
-  g->C4x[15] = (g->n*(3584*g->n-3328)+1144)/315315;
-  g->C4x[16] = (1024*g->n-208)/105105;
-  g->C4x[17] = -136/(real)(63063);
-  g->C4x[18] = (832-2560*g->n)/405405;
-  g->C4x[19] = -128/(real)(135135);
-  g->C4x[20] = 128/(real)(99099);
+  static const real coeff[] = {
+    /* C4[0], coeff of eps^5, polynomial in n of order 0 */
+    97, 15015,
+    /* C4[0], coeff of eps^4, polynomial in n of order 1 */
+    1088, 156, 45045,
+    /* C4[0], coeff of eps^3, polynomial in n of order 2 */
+    -224, -4784, 1573, 45045,
+    /* C4[0], coeff of eps^2, polynomial in n of order 3 */
+    -10656, 14144, -4576, -858, 45045,
+    /* C4[0], coeff of eps^1, polynomial in n of order 4 */
+    64, 624, -4576, 6864, -3003, 15015,
+    /* C4[0], coeff of eps^0, polynomial in n of order 5 */
+    100, 208, 572, 3432, -12012, 30030, 45045,
+    /* C4[1], coeff of eps^5, polynomial in n of order 0 */
+    1, 9009,
+    /* C4[1], coeff of eps^4, polynomial in n of order 1 */
+    -2944, 468, 135135,
+    /* C4[1], coeff of eps^3, polynomial in n of order 2 */
+    5792, 1040, -1287, 135135,
+    /* C4[1], coeff of eps^2, polynomial in n of order 3 */
+    5952, -11648, 9152, -2574, 135135,
+    /* C4[1], coeff of eps^1, polynomial in n of order 4 */
+    -64, -624, 4576, -6864, 3003, 135135,
+    /* C4[2], coeff of eps^5, polynomial in n of order 0 */
+    8, 10725,
+    /* C4[2], coeff of eps^4, polynomial in n of order 1 */
+    1856, -936, 225225,
+    /* C4[2], coeff of eps^3, polynomial in n of order 2 */
+    -8448, 4992, -1144, 225225,
+    /* C4[2], coeff of eps^2, polynomial in n of order 3 */
+    -1440, 4160, -4576, 1716, 225225,
+    /* C4[3], coeff of eps^5, polynomial in n of order 0 */
+    -136, 63063,
+    /* C4[3], coeff of eps^4, polynomial in n of order 1 */
+    1024, -208, 105105,
+    /* C4[3], coeff of eps^3, polynomial in n of order 2 */
+    3584, -3328, 1144, 315315,
+    /* C4[4], coeff of eps^5, polynomial in n of order 0 */
+    -128, 135135,
+    /* C4[4], coeff of eps^4, polynomial in n of order 1 */
+    -2560, 832, 405405,
+    /* C4[5], coeff of eps^5, polynomial in n of order 0 */
+    128, 99099,
+  };
+  int o = 0, k = 0, l, j;
+  for (l = 0; l < nC4; ++l) {        /* l is index of C4[l] */
+    for (j = nC4 - 1; j >= l; --j) { /* coeff of eps^j */
+      int m = nC4 - j - 1;           /* order of polynomial in n */
+      g->C4x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
+      o += m + 2;
+    }
+  }
 }
 
 int transit(real lon1, real lon2) {
@@ -1594,7 +1692,7 @@ void geod_polygon_addedge(const struct geod_geodesic* g,
                           real azi, real s) {
   if (p->num) {                 /* Do nothing is num is zero */
     real lat, lon, S12;
-    geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_NOWRAP, s,
+    geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s,
                    &lat, &lon, 0,
                    0, 0, 0, 0, p->polyline ? 0 : &S12);
     accadd(p->P, s);
@@ -1731,7 +1829,7 @@ unsigned geod_polygon_testedge(const struct geod_geodesic* g,
   crossings = p->crossings;
   {
     real lat, lon, s12, S12;
-    geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_NOWRAP, s,
+    geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s,
                    &lat, &lon, 0,
                    0, 0, 0, 0, &S12);
     tempsum += S12;