Issue #490 update from geodesic routines from GeographicLib 1.47.

    
Improve accuracy of area calculation (fixing a flaw introduced in
version 1.46).  Changed files geodesic.[ch3], geodtest.c, geod.1.

5 files changed, 52 insertions, 28 deletions

diff --git a/man/man1/geod.1 b/man/man1/geod.1
index 0a847f86..00400a5c 100644
--- a/man/man1/geod.1
+++ b/man/man1/geod.1
@@ -217,7 +217,7 @@ C. F. F. Karney, \fIAlgorithms for Geodesics\fR,
 .br
 J. Geodesy \fB87\fR, 43-55 (2013);
 .br
-DOI: http://dx.doi.org/10.1007/s00190-012-0578-z
+DOI: https://doi.org/10.1007/s00190-012-0578-z
 .br
 http://geographiclib.sf.net/geod-addenda.html
 .PP
diff --git a/man/man3/geodesic.3 b/man/man3/geodesic.3
index 085df834..78f6c0f3 100644
--- a/man/man3/geodesic.3
+++ b/man/man3/geodesic.3
@@ -101,7 +101,7 @@ C. F. F. Karney, \fIAlgorithms for Geodesics\fR,
 .br
 J. Geodesy \fB87\fR, 43-55 (2013);
 .br
-DOI: http://dx.doi.org/10.1007/s00190-012-0578-z
+DOI: https://doi.org/10.1007/s00190-012-0578-z
 .br
 http://geographiclib.sf.net/geod-addenda.html
 .PP
diff --git a/src/geodesic.c b/src/geodesic.c
index ad4a2686..fd964226 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -13,12 +13,12 @@
  *   C. F. F. Karney,
  *   Algorithms for geodesics,
  *   J. Geodesy <b>87</b>, 43--55 (2013);
- *   https://dx.doi.org/10.1007/s00190-012-0578-z
+ *   https://doi.org/10.1007/s00190-012-0578-z
  *   Addenda: http://geographiclib.sourceforge.net/geod-addenda.html
  *
  * See the comments in geodesic.h for documentation.
  *
- * Copyright (c) Charles Karney (2012-2016) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2017) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * http://geographiclib.sourceforge.net/
  */
@@ -274,7 +274,7 @@ static real Lambda12(const struct geod_geodesic* g,
                      real* pssig1, real* pcsig1,
                      real* pssig2, real* pcsig2,
                      real* peps,
-                     real* psomg12, real* pcomg12,
+                     real* pgomg12,
                      boolx diffp, real* pdlam12,
                      /* Scratch area of the right size */
                      real Ca[]);
@@ -786,7 +786,7 @@ static real geod_geninverse_int(const struct geod_geodesic* g,
 
     /* sig12 = sig2 - sig1 */
     sig12 = atan2(maxx((real)(0), csig1 * ssig2 - ssig1 * csig2),
-                  csig1 * csig2 + ssig1 * ssig2);
+                                  csig1 * csig2 + ssig1 * ssig2);
     Lengths(g, g->n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
             cbet1, cbet2, &s12x, &m12x, 0,
             outmask & GEOD_GEODESICSCALE ? &M12 : 0,
@@ -858,7 +858,7 @@ static real geod_geninverse_int(const struct geod_geodesic* g,
        * value of alp1 is then further from the solution) or if the new
        * estimate of alp1 lies outside (0,pi); in this case, the new starting
        * guess is taken to be (alp1a + alp1b) / 2. */
-      real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0;
+      real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0, domg12 = 0;
       unsigned numit = 0;
       /* Bracketing range */
       real salp1a = tiny, calp1a = 1, salp1b = tiny, calp1b = -1;
@@ -870,7 +870,7 @@ static real geod_geninverse_int(const struct geod_geodesic* g,
           v = Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
                         slam12, clam12,
                         &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2,
-                        &eps, &somg12, &comg12, numit < maxit1, &dv, Ca);
+                        &eps, &domg12, 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 : 1) * tol0)) break;
@@ -918,6 +918,12 @@ static real geod_geninverse_int(const struct geod_geodesic* g,
       m12x *= g->b;
       s12x *= g->b;
       a12 = sig12 / degree;
+      if (outmask & GEOD_AREA) {
+        /* omg12 = lam12 - domg12 */
+        real sdomg12 = sin(domg12), cdomg12 = cos(domg12);
+        somg12 = slam12 * cdomg12 - clam12 * sdomg12;
+        comg12 = clam12 * cdomg12 + slam12 * sdomg12;
+      }
     }
   }
 
@@ -953,11 +959,8 @@ static real geod_geninverse_int(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 && somg12 > 1) {
+      somg12 = sin(omg12); comg12 = cos(omg12);
     }
 
     if (!meridian &&
@@ -1385,15 +1388,15 @@ real Lambda12(const struct geod_geodesic* g,
               real* pssig1, real* pcsig1,
               real* pssig2, real* pcsig2,
               real* peps,
-              real* psomg12, real* pcomg12,
+              real* pdomg12,
               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,
-    somg12 = 0, comg12 = 0, dlam12 = 0;
+    domg12 = 0, dlam12 = 0;
   real salp0, calp0;
-  real somg1, comg1, somg2, comg2, lam12;
+  real somg1, comg1, somg2, comg2, somg12, comg12, lam12;
   real B312, eta, k2;
 
   if (sbet1 == 0 && calp1 == 0)
@@ -1436,7 +1439,7 @@ real Lambda12(const struct geod_geodesic* g,
 
   /* sig12 = sig2 - sig1, limit to [0, pi] */
   sig12 = atan2(maxx((real)(0), csig1 * ssig2 - ssig1 * csig2),
-                csig1 * csig2 + ssig1 * ssig2);
+                                csig1 * csig2 + ssig1 * ssig2);
 
   /* omg12 = omg2 - omg1, limit to [0, pi] */
   somg12 = maxx((real)(0), comg1 * somg2 - somg1 * comg2);
@@ -1449,7 +1452,8 @@ real Lambda12(const struct geod_geodesic* g,
   C3f(g, eps, Ca);
   B312 = (SinCosSeries(TRUE, ssig2, csig2, Ca, nC3-1) -
           SinCosSeries(TRUE, ssig1, csig1, Ca, nC3-1));
-  lam12 = eta - g->f * A3f(g, eps) * salp0 * (sig12 + B312);
+  domg12 = -g->f * A3f(g, eps) * salp0 * (sig12 + B312);
+  lam12 = eta + domg12;
 
   if (diffp) {
     if (calp2 == 0)
@@ -1469,8 +1473,7 @@ real Lambda12(const struct geod_geodesic* g,
   *pssig2 = ssig2;
   *pcsig2 = csig2;
   *peps = eps;
-  *psomg12 = somg12;
-  *pcomg12 = comg12;
+  *pdomg12 = domg12;
   if (diffp)
     *pdlam12 = dlam12;
 
@@ -1484,7 +1487,7 @@ real A3f(const struct geod_geodesic* g, real eps) {
 
 void C3f(const struct geod_geodesic* g, real eps, real c[]) {
   /* Evaluate C3 coeffs
-   * Elements c[1] through c[nC3 - 1] are set */
+   * Elements c[1] thru c[nC3 - 1] are set */
   real mult = 1;
   int o = 0, l;
   for (l = 1; l < nC3; ++l) {   /* l is index of C3[l] */
@@ -1497,7 +1500,7 @@ void C3f(const struct geod_geodesic* g, real eps, real c[]) {
 
 void C4f(const struct geod_geodesic* g, real eps, real c[]) {
   /* Evaluate C4 coeffs
-   * Elements c[0] through c[nC4 - 1] are set */
+   * Elements c[0] thru c[nC4 - 1] are set */
   real mult = 1;
   int o = 0, l;
   for (l = 0; l < nC4; ++l) {   /* l is index of C4[l] */
diff --git a/src/geodesic.h b/src/geodesic.h
index ff10c906..cd8989db 100644
--- a/src/geodesic.h
+++ b/src/geodesic.h
@@ -4,10 +4,10 @@
  *
  * This an implementation in C of the geodesic algorithms described in
  * - C. F. F. Karney,
- *   <a href="https://dx.doi.org/10.1007/s00190-012-0578-z">
+ *   <a href="https://doi.org/10.1007/s00190-012-0578-z">
  *   Algorithms for geodesics</a>,
  *   J. Geodesy <b>87</b>, 43--55 (2013);
- *   DOI: <a href="https://dx.doi.org/10.1007/s00190-012-0578-z">
+ *   DOI: <a href="https://doi.org/10.1007/s00190-012-0578-z">
  *   10.1007/s00190-012-0578-z</a>;
  *   addenda: <a href="http://geographiclib.sourceforge.net/geod-addenda.html">
  *   geod-addenda.html</a>.
@@ -112,7 +112,7 @@
  * http://geographiclib.sourceforge.net/
  *
  * This library was distributed with
- * <a href="../index.html">GeographicLib</a> 1.46.
+ * <a href="../index.html">GeographicLib</a> 1.47.
  **********************************************************************/
 
 #if !defined(GEODESIC_H)
@@ -127,12 +127,12 @@
  * The minor version of the geodesic library.  (This tracks the version of
  * GeographicLib.)
  **********************************************************************/
-#define GEODESIC_VERSION_MINOR 46
+#define GEODESIC_VERSION_MINOR 47
 /**
  * The patch level of the geodesic library.  (This tracks the version of
  * GeographicLib.)
  **********************************************************************/
-#define GEODESIC_VERSION_PATCH 1
+#define GEODESIC_VERSION_PATCH 0
 
 /**
  * Pack the version components into a single integer.  Users should not rely on
diff --git a/src/geodtest.c b/src/geodtest.c
index 990a213e..c94683f6 100644
--- a/src/geodtest.c
+++ b/src/geodtest.c
@@ -4,7 +4,7 @@
  *
  * Run these tests by configuring with cmake and running "make test".
  *
- * Copyright (c) Charles Karney (2015-2016) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2015-2017) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * http://geographiclib.sourceforge.net/
  **********************************************************************/
@@ -630,6 +630,26 @@ void polylength(const struct geod_geodesic* g, double points[][2], int N,
   geod_polygon_compute(g, &p, 0, 1, 0, perimeter);
 }
 
+int GeodSolve74() {
+  /* Check fix for inaccurate areas, bug introduced in v1.46, fixed
+     2015-10-16. */
+  double a12, s12, azi1, azi2, m12, M12, M21, S12;
+  struct geod_geodesic g;
+  int result = 0;
+  geod_init(&g, wgs84_a, wgs84_f);
+  a12 = geod_geninverse(&g, 54.1589, 15.3872, 54.1591, 15.3877,
+                        &s12, &azi1, &azi2, &m12, &M12, &M21, &S12);
+  result += assertEquals(azi1, 55.723110355, 5e-9);
+  result += assertEquals(azi2, 55.723515675, 5e-9);
+  result += assertEquals(s12,  39.527686385, 5e-9);
+  result += assertEquals(a12,   0.000355495, 5e-9);
+  result += assertEquals(m12,  39.527686385, 5e-9);
+  result += assertEquals(M12,   0.999999995, 5e-9);
+  result += assertEquals(M21,   0.999999995, 5e-9);
+  result += assertEquals(S12, 286698586.30197, 5e-4);
+  return result;
+}
+
 int Planimeter0() {
   /* Check fix for pole-encircling bug found 2011-03-16 */
   double pa[4][2] = {{89, 0}, {89, 90}, {89, 180}, {89, 270}};
@@ -757,6 +777,7 @@ int main() {
   if ((i = GeodSolve67())) {++n; printf("GeodSolve67 fail: %d\n", i);}
   if ((i = GeodSolve71())) {++n; printf("GeodSolve71 fail: %d\n", i);}
   if ((i = GeodSolve73())) {++n; printf("GeodSolve73 fail: %d\n", i);}
+  if ((i = GeodSolve74())) {++n; printf("GeodSolve74 fail: %d\n", i);}
   if ((i = Planimeter0())) {++n; printf("Planimeter0 fail: %d\n", i);}
   if ((i = Planimeter5())) {++n; printf("Planimeter5 fail: %d\n", i);}
   if ((i = Planimeter6())) {++n; printf("Planimeter6 fail: %d\n", i);}