Merge pull request #502 from cffk/geod-1.48

Merge is geodesic routines from GeographicLib 1.48.  Changes:

 * http://geographiclib.sf.net -> http://geographiclib.sourceforge.io
 * backport fixes for warnings messages from some compilers
 * change default range for longitude and azimuth to (-180d, 180d] (instead of [-180d, 180d))

plus fix compiler error so geodesic version is now 1.48.1

3 files changed, 43 insertions, 42 deletions

diff --git a/src/geodesic.c b/src/geodesic.c
index d93268c6..aeb82c71 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -14,13 +14,13 @@
  *   Algorithms for geodesics,
  *   J. Geodesy <b>87</b>, 43--55 (2013);
  *   https://doi.org/10.1007/s00190-012-0578-z
- *   Addenda: http://geographiclib.sourceforge.net/geod-addenda.html
+ *   Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
  *
  * See the comments in geodesic.h for documentation.
  *
  * Copyright (c) Charles Karney (2012-2017) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
- * http://geographiclib.sourceforge.net/
+ * https://geographiclib.sourceforge.io/
  */
 
 #include "geodesic.h"
@@ -86,8 +86,8 @@ static void Init() {
     xthresh = 1000 * tol2;
     degree = pi/180;
     {
-        double minus1 = -1.0;
-        NaN = sqrt(minus1);
+      real minus1 = -1;
+      NaN = sqrt(minus1);
     }
     init = 1;
   }
@@ -171,21 +171,21 @@ static void norm2(real* sinx, real* cosx) {
 
 static real AngNormalize(real x) {
   x = fmod(x, (real)(360));
-  return x < -180 ? x + 360 : (x < 180 ? x : x - 360);
+  return x <= -180 ? x + 360 : (x <= 180 ? x : x - 360);
 }
 
 static real LatFix(real x)
 { return fabs(x) > 90 ? NaN : x; }
 
 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
+  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
    * addition of t takes the result outside the range (-180,180] is d = 180
-   * and t < 0.  The case, d = -180 + eps, t = eps, can't happen, since
+   * 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 sumx(d == 180 && t < 0 ? -180 : d, -t, e);
+  return sumx(d == 180 && t > 0 ? -180 : d, t, e);
 }
 
 static real AngRound(real x) {
@@ -210,11 +210,12 @@ static void sincosdx(real x, real* sinx, real* cosx) {
   /* Possibly could call the gnu extension sincos */
   s = sin(r); c = cos(r);
   switch ((unsigned)q & 3U) {
-  case 0U: *sinx =     s; *cosx =     c; break;
-  case 1U: *sinx =     c; *cosx = 0 - s; break;
-  case 2U: *sinx = 0 - s; *cosx = 0 - c; break;
-  default: *sinx = 0 - c; *cosx =     s; break; /* case 3U */
+  case 0U: *sinx =  s; *cosx =  c; break;
+  case 1U: *sinx =  c; *cosx = -s; break;
+  case 2U: *sinx = -s; *cosx = -c; break;
+  default: *sinx = -c; *cosx =  s; break; /* case 3U */
   }
+  if (x != 0) { *sinx += (real)(0); *cosx += (real)(0); }
 }
 
 static real atan2dx(real y, real x) {
@@ -233,7 +234,7 @@ static real atan2dx(real y, real x) {
      *
      *   case 0: ang = 0 + ang; break;
      */
-  case 1: ang = (y > 0 ? 180 : -180) - ang; break;
+  case 1: ang = (y >= 0 ? 180 : -180) - ang; break;
   case 2: ang =  90 - ang; break;
   case 3: ang = -90 + ang; break;
   }
@@ -1752,8 +1753,8 @@ int transit(real lon1, real lon2) {
   lon1 = AngNormalize(lon1);
   lon2 = AngNormalize(lon2);
   lon12 = AngDiff(lon1, lon2, 0);
-  return lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 :
-    (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0);
+  return lon1 <= 0 && lon2 > 0 && lon12 > 0 ? 1 :
+    (lon2 <= 0 && lon1 > 0 && lon12 < 0 ? -1 : 0);
 }
 
 int transitdirect(real lon1, real lon2) {
@@ -1816,7 +1817,7 @@ void geod_polygon_addpoint(const struct geod_geodesic* g,
     p->lat0 = p->lat = lat;
     p->lon0 = p->lon = lon;
   } else {
-    real s12, S12;
+    real s12, S12 = 0;       /* Initialize S12 to stop Visual Studio warning */
     geod_geninverse(g, p->lat, p->lon, lat, lon,
                     &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12);
     accadd(p->P, s12);
@@ -1833,7 +1834,7 @@ void geod_polygon_addedge(const struct geod_geodesic* g,
                           struct geod_polygon* p,
                           real azi, real s) {
   if (p->num) {                 /* Do nothing is num is zero */
-    real lat, lon, S12;
+    real lat, lon, S12 = 0;  /* Initialize S12 to stop Visual Studio warning */
     geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s,
                    &lat, &lon, 0,
                    0, 0, 0, 0, p->polyline ? 0 : &S12);
@@ -1908,7 +1909,7 @@ unsigned geod_polygon_testpoint(const struct geod_geodesic* g,
   tempsum = p->polyline ? 0 : p->A[0];
   crossings = p->crossings;
   for (i = 0; i < (p->polyline ? 1 : 2); ++i) {
-    real s12, S12;
+    real s12, S12 = 0;       /* Initialize S12 to stop Visual Studio warning */
     geod_geninverse(g,
                     i == 0 ? p->lat  : lat, i == 0 ? p->lon  : lon,
                     i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon,
diff --git a/src/geodesic.h b/src/geodesic.h
index cd8989db..f3cb3009 100644
--- a/src/geodesic.h
+++ b/src/geodesic.h
@@ -9,7 +9,7 @@
  *   J. Geodesy <b>87</b>, 43--55 (2013);
  *   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">
+ *   addenda: <a href="https://geographiclib.sourceforge.io/geod-addenda.html">
  *   geod-addenda.html</a>.
  * .
  * The principal advantages of these algorithms over previous ones (e.g.,
@@ -96,7 +96,7 @@
  *   [\e d, \e d], for arbitrary \e d.
  *
  * These routines are a simple transcription of the corresponding C++ classes
- * in <a href="http://geographiclib.sourceforge.net"> GeographicLib</a>.  The
+ * in <a href="https://geographiclib.sourceforge.io"> GeographicLib</a>.  The
  * "class data" is represented by the structs geod_geodesic, geod_geodesicline,
  * geod_polygon and pointers to these objects are passed as initial arguments
  * to the member functions.  Most of the internal comments have been retained.
@@ -107,12 +107,12 @@
  * twice about restructuring the internals of the C code since this may make
  * porting fixes from the C++ code more difficult.
  *
- * 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/
+ * https://geographiclib.sourceforge.io/
  *
  * This library was distributed with
- * <a href="../index.html">GeographicLib</a> 1.47.
+ * <a href="../index.html">GeographicLib</a> 1.48.
  **********************************************************************/
 
 #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 47
+#define GEODESIC_VERSION_MINOR 48
 /**
  * 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
@@ -242,7 +242,7 @@ extern "C" {
    *
    * \e g must have been initialized with a call to geod_init().  \e lat1
    * should be in the range [&minus;90&deg;, 90&deg;].  The values of \e lon2
-   * and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;).  Any of
+   * and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;].  Any of
    * the "return" arguments \e plat2, etc., may be replaced by 0, if you do not
    * need some quantities computed.
    *
@@ -327,7 +327,7 @@ extern "C" {
    *
    * \e g must have been initialized with a call to geod_init().  \e lat1 and
    * \e lat2 should be in the range [&minus;90&deg;, 90&deg;].  The values of
-   * \e azi1 and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;).
+   * \e azi1 and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;].
    * Any of the "return" arguments, \e ps12, etc., may be replaced by 0, if you
    * do not need some quantities computed.
    *
@@ -525,7 +525,7 @@ extern "C" {
    *
    * \e l must have been initialized with a call, e.g., to geod_lineinit(),
    * with \e caps |= GEOD_DISTANCE_IN.  The values of \e lon2 and \e azi2
-   * returned are in the range [&minus;180&deg;, 180&deg;).  Any of the
+   * returned are in the range [&minus;180&deg;, 180&deg;].  Any of the
    * "return" arguments \e plat2, etc., may be replaced by 0, if you do not
    * need some quantities computed.
    *
@@ -594,7 +594,7 @@ extern "C" {
    *
    * \e l must have been initialized with a call to geod_lineinit() with \e
    * caps |= GEOD_DISTANCE_IN.  The value \e azi2 returned is in the range
-   * [&minus;180&deg;, 180&deg;).  Any of the "return" arguments \e plat2,
+   * [&minus;180&deg;, 180&deg;].  Any of the "return" arguments \e plat2,
    * etc., may be replaced by 0, if you do not need some quantities
    * computed.  Requesting a value which \e l is not capable of computing
    * is not an error; the corresponding argument will not be altered.
diff --git a/src/geodtest.c b/src/geodtest.c
index 7b36f64b..5ca741b1 100644
--- a/src/geodtest.c
+++ b/src/geodtest.c
@@ -6,7 +6,7 @@
  *
  * Copyright (c) Charles Karney (2015-2017) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
- * http://geographiclib.sourceforge.net/
+ * https://geographiclib.sourceforge.io/
  **********************************************************************/
 
 /** @cond SKIP */
@@ -254,7 +254,7 @@ static int GeodSolve5() {
   result += assertEquals(lat2, 90, 0.5e-5);
   if (lon2 < 0) {
     result += assertEquals(lon2, -150, 0.5e-5);
-    result += assertEquals(azi2, -180, 0.5e-5);
+    result += assertEquals(fabs(azi2), 180, 0.5e-5);
   } else {
     result += assertEquals(lon2, 30, 0.5e-5);
     result += assertEquals(azi2, 0, 0.5e-5);
@@ -340,7 +340,7 @@ static int GeodSolve14() {
   struct geod_geodesic g;
   int result = 0;
   {
-    double minus1 = -1.0;
+    double minus1 = -1;
     nan = sqrt(minus1);
   }
   geod_init(&g, wgs84_a, wgs84_f);
@@ -434,11 +434,11 @@ static int GeodSolve33() {
   result += assertEquals(s12, 19980862, 0.5);
   geod_inverse(&g, 0, 0, 0, 180, &s12, &azi1, &azi2);
   result += assertEquals(azi1, 0.00000, 0.5e-5);
-  result += assertEquals(azi2, -180.00000, 0.5e-5);
+  result += assertEquals(fabs(azi2), 180.00000, 0.5e-5);
   result += assertEquals(s12, 20003931, 0.5);
   geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2);
   result += assertEquals(azi1, 0.00000, 0.5e-5);
-  result += assertEquals(azi2, -180.00000, 0.5e-5);
+  result += assertEquals(fabs(azi2), 180.00000, 0.5e-5);
   result += assertEquals(s12, 19893357, 0.5);
   geod_init(&g, 6.4e6, 0);
   geod_inverse(&g, 0, 0, 0, 179, &s12, &azi1, &azi2);
@@ -447,11 +447,11 @@ static int GeodSolve33() {
   result += assertEquals(s12, 19994492, 0.5);
   geod_inverse(&g, 0, 0, 0, 180, &s12, &azi1, &azi2);
   result += assertEquals(azi1, 0.00000, 0.5e-5);
-  result += assertEquals(azi2, -180.00000, 0.5e-5);
+  result += assertEquals(fabs(azi2), 180.00000, 0.5e-5);
   result += assertEquals(s12, 20106193, 0.5);
   geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2);
   result += assertEquals(azi1, 0.00000, 0.5e-5);
-  result += assertEquals(azi2, -180.00000, 0.5e-5);
+  result += assertEquals(fabs(azi2), 180.00000, 0.5e-5);
   result += assertEquals(s12, 19994492, 0.5);
   geod_init(&g, 6.4e6, -1/300.0);
   geod_inverse(&g, 0, 0, 0, 179, &s12, &azi1, &azi2);
@@ -468,7 +468,7 @@ static int GeodSolve33() {
   result += assertEquals(s12, 20082617, 0.5);
   geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2);
   result += assertEquals(azi1, 0.00000, 0.5e-5);
-  result += assertEquals(azi2, -180.00000, 0.5e-5);
+  result += assertEquals(fabs(azi2), 180.00000, 0.5e-5);
   result += assertEquals(s12, 20027270, 0.5);
 
   return result;
@@ -481,7 +481,7 @@ static int GeodSolve55() {
   struct geod_geodesic g;
   int result = 0;
   {
-    double minus1 = -1.0;
+    double minus1 = -1;
     nan = sqrt(minus1);
   }
   geod_init(&g, wgs84_a, wgs84_f);
@@ -521,12 +521,12 @@ static int GeodSolve61() {
                  &lat2, &lon2, &azi2, 0, 0, 0, 0, 0);
   result += assertEquals(lat2, 45.30632, 0.5e-5);
   result += assertEquals(lon2, -180, 0.5e-5);
-  result += assertEquals(azi2, -180, 0.5e-5);
+  result += assertEquals(fabs(azi2), 180, 0.5e-5);
   geod_inverseline(&l, &g, 45, 0, 80, -0.000000000000000003, 0);
   geod_genposition(&l, flags, 1e7, &lat2, &lon2, &azi2, 0, 0, 0, 0, 0);
   result += assertEquals(lat2, 45.30632, 0.5e-5);
   result += assertEquals(lon2, -180, 0.5e-5);
-  result += assertEquals(azi2, -180, 0.5e-5);
+  result += assertEquals(fabs(azi2), 180, 0.5e-5);
   return result;
 }
 
@@ -545,7 +545,7 @@ static int GeodSolve65() {
                          &lat2, &lon2, &azi2, &s12, &m12, &M12, &M21, &S12);
   result += assertEquals(lat2, -60.23169, 0.5e-5);
   result += assertEquals(lon2, -0.00000, 0.5e-5);
-  result += assertEquals(azi2, -180.00000, 0.5e-5);
+  result += assertEquals(fabs(azi2), 180.00000, 0.5e-5);
   result += assertEquals(s12, 10000000, 0.5);
   result += assertEquals(a12, 90.06544, 0.5e-5);
   result += assertEquals(m12, 6363636, 0.5);