GeographicLib / News: Recent posts

The updates for version 1.46 introduced a small additional error in the
computation of areas. In a test case of a polygon of area 300000 km^2
with 68000 vertices representing Poland, the error is about 0.5 m^2.
The following patch corrects this problem and reduces the error to about
0.001 m^2 (comparable to the accuracy of version 1.45). This patch will
be incorporated in version 1.47

diff --git a/include/GeographicLib/Geodesic.hpp b/include/GeographicLib/Geodesic.hpp
index f932379..07390fb 100644
--- a/include/GeographicLib/Geodesic.hpp
+++ b/include/GeographicLib/Geodesic.hpp
@@ -227,7 +227,7 @@ namespace GeographicLib {
                   real salp1, real calp1, real slam120, real clam120,
                   real& salp2, real& calp2, real& sig12,
                   real& ssig1, real& csig1, real& ssig2, real& csig2,
-                  real& eps, real& somg12, real& comg12,
+                  real& eps, real& domg12,
                   bool diffp, real& dlam12, real Ca[]) const;
     real GenInverse(real lat1, real lon1, real lat2, real lon2,
                     unsigned outmask, real& s12,
diff --git a/include/GeographicLib/GeodesicExact.hpp b/include/GeographicLib/GeodesicExact.hpp
index d80098f..c8f8b3a 100644
--- a/include/GeographicLib/GeodesicExact.hpp
+++ b/include/GeographicLib/GeodesicExact.hpp
@@ -125,7 +125,7 @@ namespace GeographicLib {
                   real& salp2, real& calp2, real& sig12,
                   real& ssig1, real& csig1, real& ssig2, real& csig2,
                   EllipticFunction& E,
-                  real& somg12, real& comg12, bool diffp, real& dlam12) const;
+                  real& domg12, bool diffp, real& dlam12) const;
     real GenInverse(real lat1, real lon1, real lat2, real lon2,
                     unsigned outmask, real& s12,
                     real& salp1, real& calp1, real& salp2, real& calp2,
diff --git a/src/Geodesic.cpp b/src/Geodesic.cpp
index cf50e82..00c003a 100644
--- a/src/Geodesic.cpp
+++ b/src/Geodesic.cpp
@@ -334,7 +334,7 @@ namespace GeographicLib {
         // guess is taken to be (alp1a + alp1b) / 2.
         //
         // initial values to suppress warnings (if loop is executed 0 times)
-        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;
@@ -347,7 +347,7 @@ namespace GeographicLib {
           real v = Lambda12(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);
           // Reversed test to allow escape with NaNs
           if (tripb || !(abs(v) >= (tripn ? 8 : 1) * tol0_)) break;
           // Update bracketing values
@@ -399,6 +399,12 @@ namespace GeographicLib {
         m12x *= _b;
         s12x *= _b;
         a12 = sig12 / Math::degree();
+        if (outmask & AREA) {
+          // omg12 = lam12 - domg12
+          real sdomg12 = sin(domg12), cdomg12 = cos(domg12);
+          somg12 = slam12 * cdomg12 - clam12 * sdomg12;
+          comg12 = clam12 * cdomg12 + slam12 * sdomg12;
+        }
       }
     }

@@ -437,8 +443,7 @@ namespace GeographicLib {
       if (!meridian) {
         if (somg12 > 1) {
           somg12 = sin(omg12); comg12 = cos(omg12);
-        } else
-          Math::norm(somg12, comg12);
+        }
       }

       if (!meridian &&
@@ -810,7 +815,7 @@ namespace GeographicLib {
                                 real& sig12,
                                 real& ssig1, real& csig1,
                                 real& ssig2, real& csig2,
-                                real& eps, real& somg12, real& comg12,
+                                real& eps, real& domg12,
                                 bool diffp, real& dlam12,
                                 // Scratch area of the right size
                                 real Ca[]) const {
@@ -825,7 +830,7 @@ namespace GeographicLib {
       salp0 = salp1 * cbet1,
       calp0 = Math::hypot(calp1, salp1 * sbet1); // calp0 > 0

-    real somg1, comg1, somg2, comg2, lam12;
+    real somg1, comg1, somg2, comg2, somg12, comg12, lam12;
     // tan(bet1) = tan(sig1) * cos(alp1)
     // tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1)
     ssig1 = sbet1; somg1 = salp0 * sbet1;
@@ -871,7 +876,8 @@ namespace GeographicLib {
     C3f(eps, Ca);
     B312 = (SinCosSeries(true, ssig2, csig2, Ca, nC3_-1) -
             SinCosSeries(true, ssig1, csig1, Ca, nC3_-1));
-    lam12 = eta - _f * A3f(eps) * salp0 * (sig12 + B312);
+    domg12 = -_f * A3f(eps) * salp0 * (sig12 + B312);
+    lam12 = eta + domg12;

     if (diffp) {
       if (calp2 == 0)
diff --git a/src/GeodesicExact.cpp b/src/GeodesicExact.cpp
index 919d077..fead268 100644
--- a/src/GeodesicExact.cpp
+++ b/src/GeodesicExact.cpp
@@ -343,7 +343,7 @@ namespace GeographicLib {
         // guess is taken to be (alp1a + alp1b) / 2.
         //
         // initial values to suppress warnings (if loop is executed 0 times)
-        real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0;
+        real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, domg12 = 0;
         unsigned numit = 0;
         // Bracketing range
         real salp1a = tiny_, calp1a = 1, salp1b = tiny_, calp1b = -1;
@@ -375,7 +375,7 @@ namespace GeographicLib {
           real v = Lambda12(sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
                             slam12, clam12,
                             salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
-                            E, somg12, comg12, numit < maxit1_, dv);
+                            E, domg12, numit < maxit1_, dv);
           // Reversed test to allow escape with NaNs
           if (tripb || !(abs(v) >= (tripn ? 8 : 1) * tol0_)) break;
           // Update bracketing values
@@ -423,6 +423,12 @@ namespace GeographicLib {
         m12x *= _b;
         s12x *= _b;
         a12 = sig12 / Math::degree();
+        if (outmask & AREA) {
+          // omg12 = lam12 - domg12
+          real sdomg12 = sin(domg12), cdomg12 = cos(domg12);
+          somg12 = slam12 * cdomg12 - clam12 * sdomg12;
+          comg12 = clam12 * cdomg12 + slam12 * sdomg12;
+        }
       }
     }

@@ -462,8 +468,7 @@ namespace GeographicLib {
       if (!meridian) {
         if (somg12 > 1) {
           somg12 = sin(omg12); comg12 = cos(omg12);
-        } else
-          Math::norm(somg12, comg12);
+        }
       }

       if (!meridian &&
@@ -822,7 +827,7 @@ namespace GeographicLib {
                                      real& ssig1, real& csig1,
                                      real& ssig2, real& csig2,
                                      EllipticFunction& E,
-                                     real& somg12, real& comg12,
+                                     real& domg12,
                                      bool diffp, real& dlam12) const
     {

@@ -836,7 +841,7 @@ namespace GeographicLib {
       salp0 = salp1 * cbet1,
       calp0 = Math::hypot(calp1, salp1 * sbet1); // calp0 > 0

-    real somg1, comg1, somg2, comg2, cchi1, cchi2, lam12;
+    real somg1, comg1, somg2, comg2, somg12, comg12, cchi1, cchi2, lam12;
     // tan(bet1) = tan(sig1) * cos(alp1)
     // tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1)
     ssig1 = sbet1; somg1 = salp0 * sbet1;
@@ -888,11 +893,12 @@ namespace GeographicLib {
     // eta = chi12 - lam120
     real eta = atan2(schi12 * clam120 - cchi12 * slam120,
                      cchi12 * clam120 + schi12 * slam120);
-
-    lam12 = eta -
-      _e2/_f1 * salp0 * E.H() / (Math::pi() / 2) *
+    real deta12 = -_e2/_f1 * salp0 * E.H() / (Math::pi() / 2) *
       (sig12 + (E.deltaH(ssig2, csig2, dn2) - E.deltaH(ssig1, csig1, dn1)));
-
+    lam12 = eta + deta12;
+    // domg12 = deta12 + chi12 - omg12
+    domg12 = deta12 + atan2(schi12 * comg12 - cchi12 * somg12,
+                            cchi12 * comg12 + schi12 * somg12);
     if (diffp) {
       if (calp2 == 0)
         dlam12 = - 2 * _f1 * dn1 / sbet1;
diff --git a/tools/tests.cmake b/tools/tests.cmake
index 034e29f..31a040c 100644
--- a/tools/tests.cmake
+++ b/tools/tests.cmake
@@ -390,6 +390,17 @@ add_test (NAME GeodSolve73 COMMAND GeodSolve
 set_tests_properties (GeodSolve73
   PROPERTIES PASS_REGULAR_EXPRESSION "81\\.04623 -170\\.00000 0\\.00000")

+# Check fix for inaccurate areas, bug introduced in v1.46, fixed
+# 2015-10-16.
+add_test (NAME GeodSolve74 COMMAND GeodSolve
+  -i -p 10 -f --input-string "54.1589 15.3872 54.1591 15.3877")
+add_test (NAME GeodSolve75 COMMAND GeodSolve
+  -i -p 10 -f --input-string "54.1589 15.3872 54.1591 15.3877" -E)
+set_tests_properties (GeodSolve74 GeodSolve75
+  PROPERTIES PASS_REGULAR_EXPRESSION
+  # Exact area is 286698586.30197
+  "54.* 15.* 55\\.72311035.* 54.* 15.* 55\\.72351567.* 39\\.52768638.* 0\\.00035549.* 39\\.52768638.* 0\\.99999999.* 0\\.99999999.* 286698586\\.(29[89]|30[0-4])")
+
 # Check fix for pole-encircling bug found 2011-03-16
 add_test (NAME Planimeter0 COMMAND Planimeter
   --input-string "89 0;89 90;89 180;89 270")

Mick Killianey provided me with the following patch to the Java package
(version 1.46) which restricts the junit dependency to testing scope.
This will be incorporated in version 1.47.

diff --git a/java/pom.xml b/java/pom.xml
index 65ab122..1a11ff1 100644
--- a/java/pom.xml
+++ b/java/pom.xml
@@ -172,6 +172,7 @@
       <groupId>junit</groupId>
       <artifactId>junit</artifactId>
       <version>4.12</version>
+      <scope>test</scope>
     </dependency>
   </dependencies>

Compiling the C library, version 1.46, triggered errors with some
compilers (because declarations where not at the beginning of the
block). The following patch fixes this (and increments the version to
1.46.1). This version will be incorporated into proj.4, version 4.9.3.

The documentation is still at

http://geographiclib.sourceforge.net/1.46/C/

diff --git a/legacy/C/geodesic.c b/legacy/C/geodesic.c
index 8d9c928..e897e89 100644
--- a/legacy/C/geodesic.c
+++ b/legacy/C/geodesic.c
@@ -187,8 +187,9 @@ static real AngDiff(real x, real y, real* e) {

 static real AngRound(real x) {
   const real z = 1/(real)(16);
+  volatile real y;
   if (x == 0) return 0;
-  volatile real y = fabs(x);
+  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;
@@ -413,8 +414,8 @@ static void geod_lineinit_int(struct geod_geodesicline* l,
 void geod_lineinit(struct geod_geodesicline* l,
                    const struct geod_geodesic* g,
                    real lat1, real lon1, real azi1, unsigned caps) {
-  azi1 = AngNormalize(azi1);
   real salp1, calp1;
+  azi1 = AngNormalize(azi1);
   /* Guard against underflow in salp0 */
   sincosdx(AngRound(azi1), &salp1, &calp1);
   geod_lineinit_int(l, g, lat1, lon1, azi1, salp1, calp1, caps);
diff --git a/legacy/C/geodesic.h b/legacy/C/geodesic.h
index 80ffac9..c3f28c7 100644
--- a/legacy/C/geodesic.h
+++ b/legacy/C/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

The PyPI package, version 1.46, had a missing file. The following
patch fixes this (and increments the version to 1.46.3). This version
is now on PyPI and version 1.46 has been removed.

The documentation is still at

http://geographiclib.sourceforge.net/1.46/python/

diff --git a/python/MANIFEST.in b/python/MANIFEST.in
index 2459fc2..9561fb1 100644
--- a/python/MANIFEST.in
+++ b/python/MANIFEST.in
@@ -1 +1 @@
-# empty!
+include README.rst
diff --git a/python/geographiclib/__init__.py b/python/geographiclib/__init__.py
index 4e3507b..d889224 100644
--- a/python/geographiclib/__init__.py
+++ b/python/geographiclib/__init__.py
@@ -1,7 +1,7 @@
 """geographiclib: geodesic routines from GeographicLib"""

-__version_info__ = (1, 46, 0)
+__version_info__ = (1, 46, 3)
 """GeographicLib version as a tuple"""

-__version__ = "1.46"
+__version__ = "1.46.3"
 """GeographicLib version as a string"""
diff --git a/python/setup.py b/python/setup.py
index 2baa5b2..9a6fa69 100644
--- a/python/setup.py
+++ b/python/setup.py
@@ -38,7 +38,7 @@ class TestCommand(Command):
                                           ]))

 name = "geographiclib"
-version = "1.46"
+version = "1.46.3"

 setup(name = name,
       version = version,
@@ -46,7 +46,7 @@ setup(name = name,
       long_description = open("README.rst").read(),
       author = "Charles Karney",
       author_email = "charles@karney.com",
-      url = "http://geographiclib.sourceforge.net/" + version + "/python",
+      url = "http://geographiclib.sourceforge.net/" + "1.46" + "/python",
       packages = ["geographiclib", "geographiclib/test"],
       data_files = [],
       license = "MIT",

Changes between 1.46 (released 2016-02-15) and 1.45 versions:

• The following BUGS have been fixed:
  • the -w flag to Planimeter(1) was being ignored;
  • in the Java package, the wrong longitude was being returned with
    direct geodesic calculation with a negative distance when starting
    point was at a pole (this bug was introduced in version 1.44);
  • in the JavaScript package, PolygonArea.TestEdge contained a
    misspelling of a variable name and other typos (problem found by
    threepointone).... read more

Changes between 1.45 (released 2015-09-30) and 1.44 versions:

• Fix BUG in solution of inverse geodesic caused by misbehavior of
  some versions of Visual Studio on Windows (fmod(-0.0, 360.0) returns
  +0.0 instead of -0.0) and Octave (sind(-0.0) returns +0.0 instead of
  -0.0). These bugs were exposed because max(-0.0, +0.0) returns -0.0
  for some languages.

• Geodesic::Inverse now correctly returns NaNs if one of the latitudes
  is a NaN.... read more

The Geodesic routines are vulnerable to the "misfeature" that

max(-0.0, +0.0)

returns -0.0 with some languages. This only case of a problem caused by
this is MATLAB function geoddistance when used in Octave. The symptom
is that

geoddistance(0,0,0,[179.5,179.6])

returns two negative distances. The problem appeared in version 1.44
(released nine days ago), is only present when geoddistance is used in
Octave (not MATLAB), and is invoked with vector arguments including
nearly antipodal pairs of points on the equator.... read more

Changes between 1.44 (released 2015-08-14) and 1.43 versions:

• Various changes to improve accuracy, e.g., by minimizing round-off
  errors:
  • Add Math::sincosd, Math::sind, Math::cosd which take their
    arguments in degrees. These functions do exact range reduction
    and thus they obey exactly the elementary properties of the
    trigonometric functions, e.g., sin 9d = cos 81d =
  • sin 123456789d.
  • Math::AngNormalize now works for any angles, instead of angles in
    the range [-540d, 540d); the function Math::AngNormalize2 is now
    deprecated.
  • This means that there is now no restriction on longitudes and
    azimuths; any values can be used.
  • Improve the accuracy of Math::atan2d.
  • DMS::Decode avoids unnecessary round-off errors; thus 7:33:36 and
    7.56 result in identical values. DMS::Encode rounds ties to even.
    These changes have also been made to DMS.js.
  • More accurate rounding in MGRS::Reverse and mgrs_inv.m; this
    change only makes a difference at sub-meter precisions.
  • With MGRS::Forward and mgrs_fwd.m, ensure that digits in lower
    precision results match those at higher precision; as a result,
    strings of trailing 9s are less likely to be generated. This
    change only makes a difference at sub-meter precisions.
  • Replace the series for A2 in the Geodesic class with one with
    smaller truncation errors.
  • Geodesic::Inverse sets s12 to zero for coincident points at pole
    (instead of returning a tiny quantity).
  • Math::LatFix returns its argument if it is in [-90d, 90d]; if not,
    it returns NaN.
  • Using Math::LatFix, routines which don't check their arguments now
    interpret a latitude outside the legal range of [-90d, 90d] as a
    NaN; such routines will return NaNs instead of finite but
    incorrect results; caution: code that (dangerously) relied on the
    "reasonable" results being returned for values of the latitude
    outside the allowed range will now malfunction.... read more

There's a bug in the CassiniSoldner::Forward and cassini_fwd.m which
causes the wrong azimuth to be return for points at a pole. The
following patch fixes these problems. This will be included in the next
release.

diff --git a/src/CassiniSoldner.cpp b/src/CassiniSoldner.cpp
index f00add3..5a0f2d7 100644
--- a/src/CassiniSoldner.cpp
+++ b/src/CassiniSoldner.cpp
@@ -45,12 +45,12 @@ namespace GeographicLib {
     real sig12, s12, azi1, azi2;
     lat = Math::AngRound(lat);
     sig12 = _earth.Inverse(lat, -abs(dlon), lat, abs(dlon), s12, azi1, azi2);
-    if (sig12 < 100 * tiny_)
+    if (sig12 < 200 * tiny_)    // Need 200 > 2 / Math::degree()
       sig12 = s12 = 0;
     sig12 *= real(0.5);
     s12 *= real(0.5);
     if (s12 == 0) {
-      real da = (azi2 - azi1)/2;
+      real da = Math::AngDiff(azi1, azi2)/2;
       if (abs(dlon) <= 90) {
         azi1 = 90 - da;
         azi2 = 90 + da;
diff --git a/matlab/geographiclib/cassini_fwd.m b/matlab/geographiclib/cassini_fwd.m
index f4bf3e4..a032054 100644
--- a/matlab/geographiclib/cassini_fwd.m
+++ b/matlab/geographiclib/cassini_fwd.m
@@ -45,13 +45,13 @@ function [x, y, azi, rk] = cassini_fwd(lat0, lon0, lat, lon, ellipsoid)
   dlon = AngDiff(AngNormalize(lon0), AngNormalize(lon)) + Z;
   [s12, azi1, azi2, ~, ~, ~, ~, sig12] = ...
       geoddistance(lat, -abs(dlon), lat, abs(dlon), ellipsoid);
-  c = sig12 < 100 * tiny;
+  c = sig12 < 200 * tiny;
   sig12(c) = 0;
   s12(c) = 0;
   sig12 = 0.5 * sig12;
   s12 = 0.5 * s12;
   c = s12 == 0;
-  da = (azi2 - azi2)/2;
+  da = AngDiff(azi1, azi2)/2;
   s = abs(dlon) <= 90;
   azi1(c & s) = 90 - da(c & s);
   azi2(c & s) = 90 + da(c & s);

There are bugs in the NormalGravity class for the case of a sphere
(constructed with omega = f = J2 = 0) which caused NaNs to be returned
for most member functions. The following patch fixes these problems.
This will be included in the next release. Thanks to htallon for
finding and helping to fix these problems:

diff --git a/src/NormalGravity.cpp b/src/NormalGravity.cpp
index b673dcd..f924d86 100644
--- a/src/NormalGravity.cpp
+++ b/src/NormalGravity.cpp
@@ -30,7 +30,7 @@ namespace GeographicLib {
       if (!(Math::isfinite(_a) && _a > 0))
         throw GeographicErr("Major radius is not positive");
       if (!(Math::isfinite(_GM) && _GM > 0))
-        throw GeographicErr("Gravitational constants is not positive");
+        throw GeographicErr("Gravitational constant is not positive");
       if (!(_omega == 0 && _f == 0 && _J2 == 0)) {
         bool flatp = _f > 0 && Math::isfinite(_f);
         if (_J2 > 0 && Math::isfinite(_J2) && flatp)
@@ -40,22 +40,23 @@ namespace GeographicLib {
         if (!(Math::isfinite(_omega) && _omega != 0))
           throw GeographicErr("Angular velocity is not non-zero");
         if (flatp)
-          _J2 = FlatteningToJ2(a, GM, omega, f);
+          _J2 = FlatteningToJ2(_a, _GM, _omega, _f);
         else
-          _f = J2ToFlattening(a, GM, omega, J2);
+          _f = J2ToFlattening(_a, _GM, _omega, _J2);
       } // else we have a sphere: omega = f = J2 = 0
       _e2 = _f * (2 - _f);
       _ep2 = _e2 / (1 - _e2);
       _q0 = qf(_ep2);
       _earth = Geocentric(_a, _f);
       _b = _a * (1 - _f);
-      _E = a * sqrt(_e2);                               // H+M, Eq 2-54
-      _U0 = _GM / _E * atan(sqrt(_ep2)) + _aomega2 / 3; // H+M, Eq 2-61
+      _E = _a * sqrt(_e2);      // H+M, Eq 2-54
+      // H+M, Eq 2-61
+      _U0 = _GM / (_E ? _E / atan(sqrt(_ep2)) : _b) + _aomega2 / 3;
       // The approximate ratio of the centrifugal acceleration (at the equator)
       // to gravity.
-      _m = _aomega2 * _b / _GM;                         // H+M, Eq 2-70
+      _m = _aomega2 * _b / _GM; // H+M, Eq 2-70
       real
-        Q = _m * sqrt(_ep2) * qpf(_ep2) / (3 * _q0),
+        Q = _m * (_q0 ? sqrt(_ep2) * qpf(_ep2) / (3 * _q0) : 1),
         G = (1 - _m - Q / 2);
       _gammae = _GM / (_a * _b) * G;       // H+M, Eq 2-73
       _gammap = _GM / (_a * _a) * (1 + Q); // H+M, Eq 2-74
@@ -179,14 +180,14 @@ namespace GeographicLib {
       invw = uE / Math::hypot(u, _E * sbet), // H+M, Eq 2-63
       ep = _E/u,
       ep2 = Math::sq(ep),
-      q = qf(ep2) / _q0,
+      q = _q0 ? qf(ep2) / _q0 : pow(_a / u, 3),
       qp = qpf(ep2) / _q0,
       // H+M, Eqs 2-62 + 6-9, but omitting last (rotational) term .
-      Vres = (_GM / _E * atan(_E / u)
+      Vres = (_GM / (_E ? _E / atan(_E / u) : u)
               + _aomega2 * q * (Math::sq(sbet) - 1/real(3)) / 2),
       // H+M, Eq 6-10
       gamu = - invw * (_GM
-                       + (_aomega2 * _E * qp
+                       + (_aomega2 * (_q0 ? _E * qp : 3 * q * u)
                           * (Math::sq(sbet) - 1/real(3)) / 2)) / Math::sq(uE),
       gamb = _aomega2 * q * sbet * cbet * invw / uE,
       t = u * invw / uE;

Changes between 1.43 (released 2015-05-23) and 1.42 versions:

• Add the Enhanced Magnetic Model 2015, emm2015. This is valid for
  2000 thru the end of 2019. This required some changes in the
  MagneticModel and MagneticCircle classes; so this model cannot be
  used with versions of GeographicLib prior to 1.43.

• Fix BLUNDER in PolarStereographic constructor introduced in version
  1.42. This affected UTMUPS conversions for UPS which could be
  incorrect by up to 0.5 km.... read more

There are two serious bugs in the implementation of the LONG_NOWRAP
option for geodesics added in version 1.39. For now you should hold off
on utilizing this feature. The fixes are been made to the source tree
and will be included in version 1.43. For a preview of this release,
visit

http://geographiclib.sourceforge.net/cgi-bin/GeodSolve

and select the "Longitude: Unroll" option.

Version 1.42 introduced a bad bug into the PolarStereographic
constructor. This affected UTMUPS conversions for UPS which could be
incorrect by up to 0.5 km. The following patch fixes this. This will
be included in version 1.43.

diff --git a/src/PolarStereographic.cpp b/src/PolarStereographic.cpp
index 3ed4b1e..482251f 100644
--- a/src/PolarStereographic.cpp
+++ b/src/PolarStereographic.cpp
@@ -17,7 +17,7 @@ namespace GeographicLib {
     : _a(a)
     , _f(f <= 1 ? f : 1/f)
     , _e2(_f * (2 - _f))
-    , _es((_f < 0 ? -1 : 0) * sqrt(abs(_e2)))
+    , _es((_f < 0 ? -1 : 1) * sqrt(abs(_e2)))
     , _e2m(1 - _e2)
     , _c( (1 - _f) * exp(Math::eatanhe(real(1), _es)) )
     , _k0(k0)

Changes between 1.42 (released 2015-04-28) and 1.41 versions:

• DMS::Decode allows a single addition or subtraction operation, e.g.,
  70W+0:0:15. This affects the GeoCoords class and the utilities
  (which use the DMS class for reading coordinates).

• Add Math::norm, Math::AngRound, Math::tand, Math::atan2d,
  Math::eatanhe, Math::taupf, Math::tauf, Math::fma and remove
  duplicated (but private) functionality from other classes.... read more

The macro GEOGRAPHIC_VERSION_MINOR is incorrectly defined for autoconf
builds of GeographicLib versions 1.40 and 1.41. The problem does not
occur if building with cmake. In addition, this macro is just a
informational one; it does not affect the operation of the library. If
you are not using this macro in your code, you can safely ignore this
issue. The following patch corrects the problem for version 1.41.... read more

Changes between 1.41 (released 2015-03-09) and 1.40 versions:

• Fix bug in Rhumb::Inverse (with exact = true) and related functions
  which causes the wrong distance to be reported if one of the end
  points is at a pole. Thanks to Thomas Murray for reporting this.

• Add International Geomagnetic Reference Field (12th generation),
  which approximates the main magnetic field of the earth for the
  period 1900-2020.... read more

There's a bug in Rhumb::Inverse (with exact = true) and related
functions which causes the wrong distance to be reported if one of the
end points is at a pole.

The following patch fixes this problem. This will be included in the
next release of GeographicLib.

The bug is, in fact, in the one of the private functions of the
TransverseMercator class. But I'm not aware of an problems with the
public interface to this class.... read more

Changes between 1.40 (released 2014-12-18) and 1.39 versions:

• Add the World Magnetic Model 2015, wmm2015. This is now the default
  magnetic model for MagneticField (replacing wmm2010 which is valid
  thru the end of 2014).

• Geodesic::Inverse didn't return NaN if one of the longitudes was a
  NaN (bug introduced in version 1.25). Fixed in the C++, Java,
  JavaScript, C, Fortran, and Python implementations of the geodesic
  routines. This bug was not present in the Matlab version.... read more

I've added the World Magnetic Model 2015, wmm2015. This can be
downloaded from the magnetic-distrib folder in the download area. No
changes are needed in GeographicLib to use this model. However, I will
shortly be releasing version 1.40 in which wmm2015 will replace wmm2010
as the default magnetic model.

The following patch fixes a blunder in OSGB::GridReference (found by kalderami)
where the wrong result was returned if the easting or northing was negative. This
will be included in the next release of GeographicLib:

diff --git a/src/OSGB.cpp b/src/OSGB.cpp
index 1b83f45..d3962c7 100644
--- a/src/OSGB.cpp
+++ b/src/OSGB.cpp
@@ -44,8 +44,8 @@ namespace GeographicLib {
                           + Utility::str(int(maxprec_)) + "]");
     char grid[2 + 2 * maxprec_];
     int
-      xh = int(floor(x)) / tile_,
-      yh = int(floor(y)) / tile_;
+      xh = int(floor(x / tile_)),
+      yh = int(floor(y / tile_));
     real
       xf = x - tile_ * xh,
       yf = y - tile_ * yh;

A bug is the treatment of NaNs by Geodesic::Inverse was introduced in
version 1.25. The symptoms is that finite bogus results are returned
when a NaN is provided as one of input arguments; for example

$ echo 0 10 20 nan | GeodSolve -i
-70.06176456 -90.00000000 9987119.260

The following patch fixes this bug. This will be included in the next
release of GeographicLib.

diff --git a/src/Geodesic.cpp b/src/Geodesic.cpp
index a428ae4..d62e518 100644
--- a/src/Geodesic.cpp
+++ b/src/Geodesic.cpp
@@ -710,7 +710,8 @@ namespace GeographicLib {
         calp1 = sbet12a - cbet2 * sbet1 * Math::sq(somg12) / (1 - comg12);
       }
     }
-    if (salp1 > 0)              // Sanity check on starting guess
+    // Sanity check on starting guess.  Backwards check allows NaN through.
+    if (!(salp1 <= 0))
       SinCosNorm(salp1, calp1);
     else {
       salp1 = 1; calp1 = 0;
diff --git a/src/GeodesicExact.cpp b/src/GeodesicExact.cpp
index 47479e7..3912a24 100644
--- a/src/GeodesicExact.cpp
+++ b/src/GeodesicExact.cpp
@@ -721,7 +721,8 @@ namespace GeographicLib {
         calp1 = sbet12a - cbet2 * sbet1 * Math::sq(somg12) / (1 - comg12);
       }
     }
-    if (salp1 > 0)              // Sanity check on starting guess
+    // Sanity check on starting guess.  Backwards check allows NaN through.
+    if (!(salp1 <= 0))
       SinCosNorm(salp1, calp1);
     else {
       salp1 = 1; calp1 = 0;

Changes between 1.39 (released 2014-11-11) and 1.38 versions:

• GeographicLib usually normalizes longitudes to the range [-180deg,
  180deg). However, when solving the direct geodesic and rhumb line
  problems, it is sometimes necessary to know how many lines the line
  encircled the earth by returning the longitude "unwrapped". So the
  following changes have been made:
  • add a LONG_NOWRAP flag to mask enums for the outmask arguments for
    Geodesic, GeodesicLine, Rhumb, and RhumbLine;
  • similar changes have been made to the Python, Javascript, and
    Java implementations of the geodesic routines;
  • for the C, Fortran, and Matlab implementations the arcmode
    argument to the routines was generalized to allow a combination of
    ARCMODE and LONG_NOWRAP bits;
  • the Maxima version now returns the longitude unwrapped.
    These changes were necessary to fix the PolygonAreaT::AddEdge (see
    the next item).... read more

There's a bug in the PolygonArea::AddEdge which causes the wrong area to
be returned if the longitude extent of the added edge exceeds 180
degrees. Here's an illustration

~~~~
Geodesic geod(20e6/Math::pi(), 0);
PolygonArea poly(geod);
poly.AddPoint(0,-135);
poly.AddEdge(90,30e6);
for (int i = 135; i >= -135; --i)
poly.AddPoint(90/10e6,i);
double perimeter, area;
unsigned n = poly.Compute(false, true, perimeter, area);
cout << n << " " << perimeter << " " << area << "\n";
~~~~... read more

A bad bug was introduced into Geohash::Reverse in version 1.37. This
caused bogus results to be returned on systems where unsigned long is
only 32 bits wide. Here is the patch which will be included in version
1.39. Thanks to Christian Csar for reporting and disgnosing this
problem.

diff --git a/src/Geohash.cpp b/src/Geohash.cpp
index 00942aa..4be5594 100644
--- a/src/Geohash.cpp
+++ b/src/Geohash.cpp
@@ -89,8 +89,8 @@ namespace GeographicLib {
     int s = 5 * (maxlen_ - len);
     ulon <<=     (s / 2);
     ulat <<= s - (s / 2);
-    lon = (unsigned long)(ulon) * loneps() - 180;
-    lat = (unsigned long)(ulat) * lateps() - 90;
+    lon = ulon * loneps() - 180;
+    lat = ulat * lateps() - 90;
   }

 } // namespace GeographicLib

There's a bad bug in the area calculation in the Matlab routine
geodreckon. Probably, most people calculate areas as part of the
solution of inverse problem using geoddistance; this might explain why
this bug has gone unnoticed for so long. Here is the fix which will be
included in the next release of GeographicLib.

diff --git a/matlab/geodreckon.m b/matlab/geodreckon.m
index bdd152a..843eb0c 100644
--- a/matlab/geodreckon.m
+++ b/matlab/geodreckon.m
@@ -225,7 +225,7 @@ function [lat2, lon2, azi2, S12, m12, M12, M21, a12_s12] = geodreckon ...

   if areap
     C4x = C4coeff(n);
-    C4a = C4f(eps, C4x);
+    C4a = C4f(epsi, C4x);
     A4 = (a^2 * e2) * calp0 .* salp0;
     B41 = SinCosSeries(false, ssig1, csig1, C4a);
     B42 = SinCosSeries(false, ssig2, csig2, C4a);