Merge pull request #562 from cffk/geod-1.49

Release candidate for geodesic library version 1.49.

4 files changed, 92 insertions, 21 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index dbfefe9c..12d11c5f 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -3,6 +3,9 @@ AM_CFLAGS = @C_WFLAGS@
 bin_PROGRAMS =	proj nad2bin geod cs2cs
 EXTRA_PROGRAMS = multistresstest test228
 
+TESTS = geodtest
+check_PROGRAMS = geodtest
+
 AM_CPPFLAGS =	-DPROJ_LIB=\"$(pkgdatadir)\" \
 		-DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@
 
@@ -19,6 +22,7 @@ nad2bin_SOURCES = nad2bin.c
 geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h
 multistresstest_SOURCES = multistresstest.c
 test228_SOURCES = test228.c
+geodtest_SOURCES = geodtest.c
 
 proj_LDADD = libproj.la
 cs2cs_LDADD = libproj.la
@@ -26,6 +30,7 @@ nad2bin_LDADD = libproj.la
 geod_LDADD = libproj.la
 multistresstest_LDADD = libproj.la @THREAD_LIB@
 test228_LDADD = libproj.la @THREAD_LIB@
+geodtest_LDADD = libproj.la
 
 lib_LTLIBRARIES = libproj.la
 
diff --git a/src/geodesic.c b/src/geodesic.c
index aeb82c71..84951d7f 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -26,6 +26,10 @@
 #include "geodesic.h"
 #include <math.h>
 
+#if !defined(HAVE_C99_MATH)
+#define HAVE_C99_MATH 0
+#endif
+
 #define GEOGRAPHICLIB_GEODESIC_ORDER 6
 #define nA1   GEOGRAPHICLIB_GEODESIC_ORDER
 #define nC1   GEOGRAPHICLIB_GEODESIC_ORDER
@@ -105,6 +109,12 @@ enum captype {
 };
 
 static real sq(real x) { return x * x; }
+#if HAVE_C99_MATH
+#define atanhx atanh
+#define copysignx copysign
+#define hypotx hypot
+#define cbrtx cbrt
+#else
 static real log1px(real x) {
   volatile real
     y = 1 + x,
@@ -133,6 +143,7 @@ static real cbrtx(real x) {
   real y = pow(fabs(x), 1/(real)(3)); /* Return the real cube root */
   return x < 0 ? -y : y;
 }
+#endif
 
 static real sumx(real u, real v, real* t) {
   volatile real s = u + v;
@@ -170,8 +181,13 @@ static void norm2(real* sinx, real* cosx) {
 }
 
 static real AngNormalize(real x) {
+#if HAVE_C99_MATH
+  x = remainder(x, (real)(360));
+  return x != -180 ? x : 180;
+#else
   x = fmod(x, (real)(360));
   return x <= -180 ? x + 360 : (x <= 180 ? x : x - 360);
+#endif
 }
 
 static real LatFix(real x)
@@ -202,9 +218,15 @@ static void sincosdx(real x, real* sinx, real* cosx) {
   /* In order to minimize round-off errors, this function exactly reduces
    * the argument to the range [-45, 45] before converting it to radians. */
   real r, s, c; int q;
+#if HAVE_C99_MATH && !defined(__GNUC__)
+  /* Disable for gcc because of bug in glibc version < 2.22, see
+   * https://sourceware.org/bugzilla/show_bug.cgi?id=17569 */
+  r = remquo(x, (real)(90), &q);
+#else
   r = fmod(x, (real)(360));
   q = (int)(floor(r / 90 + (real)(0.5)));
   r -= 90 * q;
+#endif
   /* now abs(r) <= 45 */
   r *= degree;
   /* Possibly could call the gnu extension sincos */
@@ -538,7 +560,9 @@ real geod_genposition(const struct geod_geodesicline* l,
   salp2 = l->salp0; calp2 = l->calp0 * csig2; /* No need to normalize */
 
   if (outmask & GEOD_DISTANCE)
-    s12 = flags & GEOD_ARCMODE ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12;
+    s12 = flags & GEOD_ARCMODE ?
+      l->b * ((1 + l->A1m1) * sig12 + AB1) :
+      s12_a12;
 
   if (outmask & GEOD_LONGITUDE) {
     real E = copysignx(1, l->salp0); /* east or west going? */
@@ -576,7 +600,8 @@ real geod_genposition(const struct geod_geodesicline* l,
       m12 = l->b * ((dn2 * (l->csig1 * ssig2) - l->dn1 * (l->ssig1 * csig2))
                     - l->csig1 * csig2 * J12);
     if (outmask & GEOD_GEODESICSCALE) {
-      real t = l->k2 * (ssig2 - l->ssig1) * (ssig2 + l->ssig1) / (l->dn1 + dn2);
+      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;
     }
@@ -639,7 +664,9 @@ static void geod_setarc(struct geod_geodesicline* l, real a13) {
 
 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);
+  flags & GEOD_ARCMODE ?
+    geod_setarc(l, s13_a13) :
+    geod_setdistance(l, s13_a13);
 }
 
 void geod_position(const struct geod_geodesicline* l, real s12,
@@ -1758,10 +1785,17 @@ int transit(real lon1, real lon2) {
 }
 
 int transitdirect(real lon1, real lon2) {
+#if HAVE_C99_MATH
+  lon1 = remainder(lon1, (real)(720));
+  lon2 = remainder(lon2, (real)(720));
+  return ( (lon2 >= 0 && lon2 < 360 ? 0 : 1) -
+           (lon1 >= 0 && lon1 < 360 ? 0 : 1) );
+#else
   lon1 = fmod(lon1, (real)(720));
   lon2 = fmod(lon2, (real)(720));
   return ( ((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) -
            ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1) );
+#endif
 }
 
 void accini(real s[]) {
diff --git a/src/geodesic.h b/src/geodesic.h
index f3cb3009..ab18a01f 100644
--- a/src/geodesic.h
+++ b/src/geodesic.h
@@ -112,7 +112,7 @@
  * https://geographiclib.sourceforge.io/
  *
  * This library was distributed with
- * <a href="../index.html">GeographicLib</a> 1.48.
+ * <a href="../index.html">GeographicLib</a> 1.49.
  **********************************************************************/
 
 #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 48
+#define GEODESIC_VERSION_MINOR 49
 /**
  * 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
@@ -881,16 +881,16 @@ extern "C" {
    * mask values for the \e caps argument to geod_lineinit().
    **********************************************************************/
   enum geod_mask {
-    GEOD_NONE         = 0U,                     /**< Calculate nothing */
-    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 */
+    GEOD_NONE         = 0U,                    /**< Calculate nothing */
+    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 */
   };
 
   /**
diff --git a/src/geodtest.c b/src/geodtest.c
index 5ca741b1..6899436c 100644
--- a/src/geodtest.c
+++ b/src/geodtest.c
@@ -16,7 +16,7 @@
 #include <math.h>
 
 #if defined(_MSC_VER)
-// Squelch warnings about assignment within conditional expression
+/* Squelch warnings about assignment within conditional expression */
 #  pragma warning (disable: 4706)
 #endif
 
@@ -618,8 +618,9 @@ static int GeodSolve73() {
   return result;
 }
 
-static void planimeter(const struct geod_geodesic* g, double points[][2], int N,
-                double* perimeter, double* area) {
+static void planimeter(const struct geod_geodesic* g,
+                       double points[][2], int N,
+                       double* perimeter, double* area) {
   struct geod_polygon p;
   int i;
   geod_polygon_init(&p, 0);
@@ -628,8 +629,9 @@ static void planimeter(const struct geod_geodesic* g, double points[][2], int N,
   geod_polygon_compute(g, &p, 0, 1, area, perimeter);
 }
 
-static void polylength(const struct geod_geodesic* g, double points[][2], int N,
-                double* perimeter) {
+static void polylength(const struct geod_geodesic* g,
+                       double points[][2], int N,
+                       double* perimeter) {
   struct geod_polygon p;
   int i;
   geod_polygon_init(&p, 1);
@@ -658,6 +660,34 @@ static int GeodSolve74() {
   return result;
 }
 
+static int GeodSolve76() {
+  /* The distance from Wellington and Salamanca (a classic failure of
+     Vincenty) */
+  double azi1, azi2, s12;
+  struct geod_geodesic g;
+  int result = 0;
+  geod_init(&g, wgs84_a, wgs84_f);
+  geod_inverse(&g, -(41+19/60.0), 174+49/60.0, 40+58/60.0, -(5+30/60.0),
+               &s12, &azi1, &azi2);
+  result += assertEquals(azi1, 160.39137649664, 0.5e-11);
+  result += assertEquals(azi2,  19.50042925176, 0.5e-11);
+  result += assertEquals(s12,  19960543.857179, 0.5e-6);
+  return result;
+}
+
+static int GeodSolve78() {
+  /* An example where the NGS calculator fails to converge */
+  double azi1, azi2, s12;
+  struct geod_geodesic g;
+  int result = 0;
+  geod_init(&g, wgs84_a, wgs84_f);
+  geod_inverse(&g, 27.2, 0.0, -27.1, 179.5, &s12, &azi1, &azi2);
+  result += assertEquals(azi1,  45.82468716758, 0.5e-11);
+  result += assertEquals(azi2, 134.22776532670, 0.5e-11);
+  result += assertEquals(s12,  19974354.765767, 0.5e-6);
+  return result;
+}
+
 static int Planimeter0() {
   /* Check fix for pole-encircling bug found 2011-03-16 */
   double pa[4][2] = {{89, 0}, {89, 90}, {89, 180}, {89, 270}};
@@ -786,6 +816,8 @@ int main() {
   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 = GeodSolve76())) {++n; printf("GeodSolve76 fail: %d\n", i);}
+  if ((i = GeodSolve78())) {++n; printf("GeodSolve78 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);}