Fix issue OSGeo#715. Merge 1.49.1 of geodesic library (tagged as v1.4…

…9.1-c

in GeographicLib). Details:

Workaround bugs in handling of -0.0 in fmod and sin in Visual Studio 10,
11, and 12.

Relax unrealistically strict delta for GeodSolve59 in geodtest.

src/geodesic.c

@@ -185,8 +185,18 @@ static real AngNormalize(real x) {
   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);
+  real y = fmod(x, (real)(360));
+#if defined(_MSC_VER) && _MSC_VER < 1900
+  /*
+    Before version 14 (2015), Visual Studio had problems dealing
+    with -0.0.  Specifically
+      VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0
+    sincosdx has a similar fix.
+    python 2.7 on Windows 32-bit machines has the same problem.
+  */
+  if (x == 0) y = x;
+#endif
+  return y <= -180 ? y + 360 : (y <= 180 ? y : y - 360);
 #endif
 }
 
@@ -231,6 +241,17 @@ static void sincosdx(real x, real* sinx, real* cosx) {
   r *= degree;
   /* Possibly could call the gnu extension sincos */
   s = sin(r); c = cos(r);
+#if defined(_MSC_VER) && _MSC_VER < 1900
+  /*
+    Before version 14 (2015), Visual Studio had problems dealing
+    with -0.0.  Specifically
+      VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0
+      VC 12       and 64-bit compile:  sin(-0.0)        -> +0.0
+    AngNormalize has a similar fix.
+    python 2.7 on Windows 32-bit machines has the same problem.
+  */
+  if (x == 0) s = x;
+#endif
   switch ((unsigned)q & 3U) {
   case 0U: *sinx =  s; *cosx =  c; break;
   case 1U: *sinx =  c; *cosx = -s; break;

src/geodesic.h

@@ -132,7 +132,7 @@
  * The patch level of the geodesic library.  (This tracks the version of
  * GeographicLib.)
  **********************************************************************/
-#define GEODESIC_VERSION_PATCH 0
+#define GEODESIC_VERSION_PATCH 1
 
 /**
  * Pack the version components into a single integer.  Users should not rely on

src/geodtest.c

@@ -505,7 +505,7 @@ static int GeodSolve59() {
   geod_inverse(&g, 5, 0.00000000000001, 10, 180, &s12, &azi1, &azi2);
   result += assertEquals(azi1, 0.000000000000035, 1.5e-14);
   result += assertEquals(azi2, 179.99999999999996, 1.5e-14);
-  result += assertEquals(s12, 18345191.174332713, 2.5e-9);
+  result += assertEquals(s12, 18345191.174332713, 5e-9);
   return result;
 }