diff options
| -rw-r--r-- | docs/source/geodesic.rst | 4 | ||||
| -rw-r--r-- | src/geodesic.c | 9 | ||||
| -rw-r--r-- | src/geodesic.h | 6 | ||||
| -rw-r--r-- | src/tests/geodtest.cpp | 24 |
4 files changed, 33 insertions, 10 deletions
diff --git a/docs/source/geodesic.rst b/docs/source/geodesic.rst index 11bb5e5c..6d07352c 100644 --- a/docs/source/geodesic.rst +++ b/docs/source/geodesic.rst @@ -47,7 +47,7 @@ Traditionally two geodesic problems are considered: :math:`\alpha_2`. PROJ incorporates `C library for Geodesics -<https://geographiclib.sourceforge.io/1.49/C/>`_ from `GeographicLib +<https://geographiclib.sourceforge.io/1.52/C/>`_ from `GeographicLib <https://geographiclib.sourceforge.io>`_. This library provides routines to solve the direct and inverse geodesic problems. Full double precision accuracy is maintained provided that @@ -55,7 +55,7 @@ precision accuracy is maintained provided that to the `application programming interface - <https://geographiclib.sourceforge.io/1.49/C/geodesic_8h.html>`_ + <https://geographiclib.sourceforge.io/1.52/C/geodesic_8h.html>`_ for full documentation. A brief summary of the routines is given in geodesic(3). diff --git a/src/geodesic.c b/src/geodesic.c index 53ec9ed6..d5b22ac6 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -18,7 +18,7 @@ * * See the comments in geodesic.h for documentation. * - * Copyright (c) Charles Karney (2012-2020) <charles@karney.com> and licensed + * Copyright (c) Charles Karney (2012-2021) <charles@karney.com> and licensed * under the MIT/X11 License. For more information, see * https://geographiclib.sourceforge.io/ */ @@ -164,7 +164,7 @@ static real AngRound(real x) { 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; + real r, s, c; int q = 0; r = remquo(x, (real)(90), &q); /* now abs(r) <= 45 */ r *= degree; @@ -797,7 +797,9 @@ static real geod_geninverse_int(const struct geod_geodesic* g, * not a shortest path. */ if (sig12 < 1 || m12x >= 0) { /* Need at least 2, to handle 90 0 90 180 */ - if (sig12 < 3 * tiny) + if (sig12 < 3 * tiny || + // Prevent negative s12 or m12 for short lines + (sig12 < tol0 && (s12x < 0 || m12x < 0))) sig12 = m12x = s12x = 0; m12x *= g->b; s12x *= g->b; @@ -868,7 +870,6 @@ static real geod_geninverse_int(const struct geod_geodesic* g, slam12, clam12, &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2, &eps, &domg12, numit < maxit1, &dv, Ca); - /* 2 * tol0 is approximately 1 ulp for a number in [0, pi]. */ /* Reversed test to allow escape with NaNs */ if (tripb || !(fabs(v) >= (tripn ? 8 : 1) * tol0)) break; /* Update bracketing values */ diff --git a/src/geodesic.h b/src/geodesic.h index 548fbf47..b803d2a5 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -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-2020) <charles@karney.com> and licensed + * Copyright (c) Charles Karney (2012-2021) <charles@karney.com> and licensed * under the MIT/X11 License. For more information, see * https://geographiclib.sourceforge.io/ * * This library was distributed with - * <a href="../index.html">GeographicLib</a> 1.51. + * <a href="../index.html">GeographicLib</a> 1.52. **********************************************************************/ #if !defined(GEODESIC_H) @@ -127,7 +127,7 @@ * The minor version of the geodesic library. (This tracks the version of * GeographicLib.) **********************************************************************/ -#define GEODESIC_VERSION_MINOR 51 +#define GEODESIC_VERSION_MINOR 52 /** * The patch level of the geodesic library. (This tracks the version of * GeographicLib.) diff --git a/src/tests/geodtest.cpp b/src/tests/geodtest.cpp index 1b1e7aa7..c4d8b098 100644 --- a/src/tests/geodtest.cpp +++ b/src/tests/geodtest.cpp @@ -4,7 +4,7 @@ * * Run these tests by configuring with cmake and running "make test". * - * Copyright (c) Charles Karney (2015-2019) <charles@karney.com> and licensed + * Copyright (c) Charles Karney (2015-2021) <charles@karney.com> and licensed * under the MIT/X11 License. For more information, see * https://geographiclib.sourceforge.io/ **********************************************************************/ @@ -737,6 +737,9 @@ static int GeodSolve80() { result += checkEquals(M12, 1, 1e-15); result += checkEquals(M21, 1, 1e-15); result += checkEquals(S12, 0, 1e-10); + result += 1/a12 > 0 ? 0 : 1; + result += 1/s12 > 0 ? 0 : 1; + result += 1/m12 > 0 ? 0 : 1; a12 = geod_geninverse(&g, 90, 0, 90, 180, &s12, &azi1, &azi2, &m12, &M12, &M21, &S12); @@ -806,6 +809,24 @@ static int GeodSolve84() { return result; } +static int GeodSolve92() { + /* Check fix for inaccurate hypot with python 3.[89]. Problem reported + * by agdhruv https://github.com/geopy/geopy/issues/466 ; see + * https://bugs.python.org/issue43088 */ + double azi1, azi2, s12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverse(&g, + 37.757540000000006, -122.47018, + 37.75754, -122.470177, + &s12, &azi1, &azi2); + result += checkEquals(azi1, 89.99999923, 1e-7 ); + result += checkEquals(azi2, 90.00000106, 1e-7 ); + result += checkEquals(s12, 0.264, 0.5e-3); + 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}}; @@ -1092,6 +1113,7 @@ int main() { if ((i = GeodSolve78())) {++n; printf("GeodSolve78 fail: %d\n", i);} if ((i = GeodSolve80())) {++n; printf("GeodSolve80 fail: %d\n", i);} if ((i = GeodSolve84())) {++n; printf("GeodSolve84 fail: %d\n", i);} + if ((i = GeodSolve92())) {++n; printf("GeodSolve92 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);} |