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:
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:
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";
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);
The file src/GeographicLib.pro does not include all the files in
GeographicLib. This problem is fixed with the attached patch. This fix
will be included in the next release of GeographicLib.
diff --git a/src/GeographicLib.pro b/src/GeographicLib.pro index 1f3d036..332db12 100644 --- a/src/GeographicLib.pro +++ b/src/GeographicLib.pro @@ -5,6 +5,7 @@ TEMPLATE = lib INCLUDEPATH = ../include INCLUDEDIR = $$INCLUDEPATH/GeographicLib +SOURCES += Accumulator.cpp SOURCES += AlbersEqualArea.cpp SOURCES += AzimuthalEquidistant.cpp SOURCES += CassiniSoldner.cpp @@ -15,7 +16,10 @@ SOURCES += EllipticFunction.cpp SOURCES += GeoCoords.cpp SOURCES += Geocentric.cpp SOURCES += Geodesic.cpp +SOURCES += GeodesicExact.cpp +SOURCES += GeodesicExactC4.cpp SOURCES += GeodesicLine.cpp +SOURCES += GeodesicLineExact.cpp SOURCES += Geohash.cpp SOURCES += Geoid.cpp SOURCES += Gnomonic.cpp @@ -49,7 +53,9 @@ HEADERS += $$INCLUDEDIR/EllipticFunction.hpp HEADERS += $$INCLUDEDIR/GeoCoords.hpp HEADERS += $$INCLUDEDIR/Geocentric.hpp HEADERS += $$INCLUDEDIR/Geodesic.hpp +HEADERS += $$INCLUDEDIR/GeodesicExact.hpp HEADERS += $$INCLUDEDIR/GeodesicLine.hpp +HEADERS += $$INCLUDEDIR/GeodesicLineExact.hpp HEADERS += $$INCLUDEDIR/Geohash.hpp HEADERS += $$INCLUDEDIR/Geoid.hpp HEADERS += $$INCLUDEDIR/Gnomonic.hpp
Changes between 1.38 (released 2014-10-02) and 1.37 versions:
On MacOSX, the installed package is relocatable (for cmake version
2.8.12 and later).
On Mac OSX, GeographicLib can be installed using homebrew.
In cmake builds under Windows, set the output directories so that
binaries and shared libraries are together.
Accept the minus sign as a synomym for - in DMS.{cpp,js}.
The cmake configuration file geographiclib-depends.cmake has been
renamed to geographiclib-targets.cmake.... read more
Changes between 1.37 (released 2014-08-08) and 1.36 versions:
Add support for high precision arithmetic.
INCOMPATIBLE CHANGE: the static instantiations of various classes
for the WGS84 ellipsoid have been changed to a "construct on first
use idiom". This avoids a lot of wasteful initialization before the
user's code starts. Unfortunately it means that existing source
code that relies on any of the following static variables will need
to be changed to a function call:
Changes between 1.36 (released 2014-05-13) and 1.35 versions:
After upgrading to GeographicLib version 1.34 or later, it's possible
that projects that compile and link to GeographicLib will get this
compiler error:
#error : GEOGRAPHICLIB_SHARED_LIB must be 0 or 1
This will only happen on non-Windows systems and is caused by cmake
finding the cmake configuration for a pre-1.34 version of GeographicLib.
The fix is to remove all earlier cmake configuration files with (as
root):... read more