From 2a13208980a42b645e72486109a9ae02a92238a3 Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Tue, 29 Aug 2017 15:58:00 -0400 Subject: Release candidate for geodesic library version 1.49. Only substantial changes are (1) testing the HAVE_C99_MATH flag and acting accordingly and (2) adding a couple of tests. --- src/geodesic.c | 40 +++++++++++++++++++++++++++++++++++++--- src/geodesic.h | 26 +++++++++++++------------- src/geodtest.c | 42 +++++++++++++++++++++++++++++++++++++----- 3 files changed, 87 insertions(+), 21 deletions(-) (limited to 'src') diff --git a/src/geodesic.c b/src/geodesic.c index aeb82c71..84951d7f 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -26,6 +26,10 @@ #include "geodesic.h" #include +#if !defined(HAVE_C99_MATH) +#define HAVE_C99_MATH 0 +#endif + #define GEOGRAPHICLIB_GEODESIC_ORDER 6 #define nA1 GEOGRAPHICLIB_GEODESIC_ORDER #define nC1 GEOGRAPHICLIB_GEODESIC_ORDER @@ -105,6 +109,12 @@ enum captype { }; static real sq(real x) { return x * x; } +#if HAVE_C99_MATH +#define atanhx atanh +#define copysignx copysign +#define hypotx hypot +#define cbrtx cbrt +#else static real log1px(real x) { volatile real y = 1 + x, @@ -133,6 +143,7 @@ static real cbrtx(real x) { real y = pow(fabs(x), 1/(real)(3)); /* Return the real cube root */ return x < 0 ? -y : y; } +#endif static real sumx(real u, real v, real* t) { volatile real s = u + v; @@ -170,8 +181,13 @@ static void norm2(real* sinx, real* cosx) { } static real AngNormalize(real x) { +#if HAVE_C99_MATH + 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); +#endif } static real LatFix(real x) @@ -202,9 +218,15 @@ 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; +#if HAVE_C99_MATH && !defined(__GNUC__) + /* Disable for gcc because of bug in glibc version < 2.22, see + * https://sourceware.org/bugzilla/show_bug.cgi?id=17569 */ + r = remquo(x, (real)(90), &q); +#else r = fmod(x, (real)(360)); q = (int)(floor(r / 90 + (real)(0.5))); r -= 90 * q; +#endif /* now abs(r) <= 45 */ r *= degree; /* Possibly could call the gnu extension sincos */ @@ -538,7 +560,9 @@ real geod_genposition(const struct geod_geodesicline* l, salp2 = l->salp0; calp2 = l->calp0 * csig2; /* No need to normalize */ if (outmask & GEOD_DISTANCE) - s12 = flags & GEOD_ARCMODE ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12; + s12 = flags & GEOD_ARCMODE ? + l->b * ((1 + l->A1m1) * sig12 + AB1) : + s12_a12; if (outmask & GEOD_LONGITUDE) { real E = copysignx(1, l->salp0); /* east or west going? */ @@ -576,7 +600,8 @@ real geod_genposition(const struct geod_geodesicline* l, m12 = l->b * ((dn2 * (l->csig1 * ssig2) - l->dn1 * (l->ssig1 * csig2)) - l->csig1 * csig2 * J12); if (outmask & GEOD_GEODESICSCALE) { - real t = l->k2 * (ssig2 - l->ssig1) * (ssig2 + l->ssig1) / (l->dn1 + dn2); + real t = l->k2 * (ssig2 - l->ssig1) * (ssig2 + l->ssig1) / + (l->dn1 + dn2); M12 = csig12 + (t * ssig2 - csig2 * J12) * l->ssig1 / l->dn1; M21 = csig12 - (t * l->ssig1 - l->csig1 * J12) * ssig2 / dn2; } @@ -639,7 +664,9 @@ static void geod_setarc(struct geod_geodesicline* l, real a13) { void geod_gensetdistance(struct geod_geodesicline* l, unsigned flags, real s13_a13) { - flags & GEOD_ARCMODE ? geod_setarc(l, s13_a13) : geod_setdistance(l, s13_a13); + flags & GEOD_ARCMODE ? + geod_setarc(l, s13_a13) : + geod_setdistance(l, s13_a13); } void geod_position(const struct geod_geodesicline* l, real s12, @@ -1758,10 +1785,17 @@ int transit(real lon1, real lon2) { } int transitdirect(real lon1, real lon2) { +#if HAVE_C99_MATH + lon1 = remainder(lon1, (real)(720)); + lon2 = remainder(lon2, (real)(720)); + return ( (lon2 >= 0 && lon2 < 360 ? 0 : 1) - + (lon1 >= 0 && lon1 < 360 ? 0 : 1) ); +#else lon1 = fmod(lon1, (real)(720)); lon2 = fmod(lon2, (real)(720)); return ( ((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) - ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1) ); +#endif } void accini(real s[]) { diff --git a/src/geodesic.h b/src/geodesic.h index f3cb3009..ab18a01f 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -112,7 +112,7 @@ * https://geographiclib.sourceforge.io/ * * This library was distributed with - * GeographicLib 1.48. + * GeographicLib 1.49. **********************************************************************/ #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 48 +#define GEODESIC_VERSION_MINOR 49 /** * The patch level of the geodesic library. (This tracks the version of * GeographicLib.) **********************************************************************/ -#define GEODESIC_VERSION_PATCH 1 +#define GEODESIC_VERSION_PATCH 0 /** * Pack the version components into a single integer. Users should not rely on @@ -881,16 +881,16 @@ extern "C" { * mask values for the \e caps argument to geod_lineinit(). **********************************************************************/ enum geod_mask { - GEOD_NONE = 0U, /**< Calculate nothing */ - GEOD_LATITUDE = 1U<<7 | 0U, /**< Calculate latitude */ - GEOD_LONGITUDE = 1U<<8 | 1U<<3, /**< Calculate longitude */ - GEOD_AZIMUTH = 1U<<9 | 0U, /**< Calculate azimuth */ - GEOD_DISTANCE = 1U<<10 | 1U<<0, /**< Calculate distance */ - GEOD_DISTANCE_IN = 1U<<11 | 1U<<0 | 1U<<1, /**< Allow distance as input */ - GEOD_REDUCEDLENGTH= 1U<<12 | 1U<<0 | 1U<<2, /**< Calculate reduced length */ - GEOD_GEODESICSCALE= 1U<<13 | 1U<<0 | 1U<<2, /**< Calculate geodesic scale */ - GEOD_AREA = 1U<<14 | 1U<<4, /**< Calculate reduced length */ - GEOD_ALL = 0x7F80U| 0x1FU /**< Calculate everything */ + GEOD_NONE = 0U, /**< Calculate nothing */ + GEOD_LATITUDE = 1U<<7 | 0U, /**< Calculate latitude */ + GEOD_LONGITUDE = 1U<<8 | 1U<<3, /**< Calculate longitude */ + GEOD_AZIMUTH = 1U<<9 | 0U, /**< Calculate azimuth */ + GEOD_DISTANCE = 1U<<10 | 1U<<0, /**< Calculate distance */ + GEOD_DISTANCE_IN = 1U<<11 | 1U<<0 | 1U<<1,/**< Allow distance as input */ + GEOD_REDUCEDLENGTH= 1U<<12 | 1U<<0 | 1U<<2,/**< Calculate reduced length */ + GEOD_GEODESICSCALE= 1U<<13 | 1U<<0 | 1U<<2,/**< Calculate geodesic scale */ + GEOD_AREA = 1U<<14 | 1U<<4, /**< Calculate reduced length */ + GEOD_ALL = 0x7F80U| 0x1FU /**< Calculate everything */ }; /** diff --git a/src/geodtest.c b/src/geodtest.c index 5ca741b1..6899436c 100644 --- a/src/geodtest.c +++ b/src/geodtest.c @@ -16,7 +16,7 @@ #include #if defined(_MSC_VER) -// Squelch warnings about assignment within conditional expression +/* Squelch warnings about assignment within conditional expression */ # pragma warning (disable: 4706) #endif @@ -618,8 +618,9 @@ static int GeodSolve73() { return result; } -static void planimeter(const struct geod_geodesic* g, double points[][2], int N, - double* perimeter, double* area) { +static void planimeter(const struct geod_geodesic* g, + double points[][2], int N, + double* perimeter, double* area) { struct geod_polygon p; int i; geod_polygon_init(&p, 0); @@ -628,8 +629,9 @@ static void planimeter(const struct geod_geodesic* g, double points[][2], int N, geod_polygon_compute(g, &p, 0, 1, area, perimeter); } -static void polylength(const struct geod_geodesic* g, double points[][2], int N, - double* perimeter) { +static void polylength(const struct geod_geodesic* g, + double points[][2], int N, + double* perimeter) { struct geod_polygon p; int i; geod_polygon_init(&p, 1); @@ -658,6 +660,34 @@ static int GeodSolve74() { return result; } +static int GeodSolve76() { + /* The distance from Wellington and Salamanca (a classic failure of + Vincenty) */ + double azi1, azi2, s12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverse(&g, -(41+19/60.0), 174+49/60.0, 40+58/60.0, -(5+30/60.0), + &s12, &azi1, &azi2); + result += assertEquals(azi1, 160.39137649664, 0.5e-11); + result += assertEquals(azi2, 19.50042925176, 0.5e-11); + result += assertEquals(s12, 19960543.857179, 0.5e-6); + return result; +} + +static int GeodSolve78() { + /* An example where the NGS calculator fails to converge */ + double azi1, azi2, s12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverse(&g, 27.2, 0.0, -27.1, 179.5, &s12, &azi1, &azi2); + result += assertEquals(azi1, 45.82468716758, 0.5e-11); + result += assertEquals(azi2, 134.22776532670, 0.5e-11); + result += assertEquals(s12, 19974354.765767, 0.5e-6); + 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}}; @@ -786,6 +816,8 @@ int main() { if ((i = GeodSolve71())) {++n; printf("GeodSolve71 fail: %d\n", i);} if ((i = GeodSolve73())) {++n; printf("GeodSolve73 fail: %d\n", i);} if ((i = GeodSolve74())) {++n; printf("GeodSolve74 fail: %d\n", i);} + if ((i = GeodSolve76())) {++n; printf("GeodSolve76 fail: %d\n", i);} + if ((i = GeodSolve78())) {++n; printf("GeodSolve78 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);} -- cgit v1.2.3 From 03eb6a7804019df8c90014c1fb5dd885bc08b40e Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Wed, 30 Aug 2017 09:54:41 -0400 Subject: Add geodtest to autoconf's 'make check' suite. It's already including in cmake's test suite. --- src/Makefile.am | 5 +++++ 1 file changed, 5 insertions(+) (limited to 'src') diff --git a/src/Makefile.am b/src/Makefile.am index dbfefe9c..12d11c5f 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -3,6 +3,9 @@ AM_CFLAGS = @C_WFLAGS@ bin_PROGRAMS = proj nad2bin geod cs2cs EXTRA_PROGRAMS = multistresstest test228 +TESTS = geodtest +check_PROGRAMS = geodtest + AM_CPPFLAGS = -DPROJ_LIB=\"$(pkgdatadir)\" \ -DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@ @@ -19,6 +22,7 @@ nad2bin_SOURCES = nad2bin.c geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h multistresstest_SOURCES = multistresstest.c test228_SOURCES = test228.c +geodtest_SOURCES = geodtest.c proj_LDADD = libproj.la cs2cs_LDADD = libproj.la @@ -26,6 +30,7 @@ nad2bin_LDADD = libproj.la geod_LDADD = libproj.la multistresstest_LDADD = libproj.la @THREAD_LIB@ test228_LDADD = libproj.la @THREAD_LIB@ +geodtest_LDADD = libproj.la lib_LTLIBRARIES = libproj.la -- cgit v1.2.3 From f658c3c822e65ec907e44f14c937aef7824d404c Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Wed, 30 Aug 2017 10:11:33 -0400 Subject: Rejig how geodtest is specified in automake to make coverage stats accurate?? --- src/Makefile.am | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) (limited to 'src') diff --git a/src/Makefile.am b/src/Makefile.am index 12d11c5f..d215e678 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -22,7 +22,7 @@ nad2bin_SOURCES = nad2bin.c geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h multistresstest_SOURCES = multistresstest.c test228_SOURCES = test228.c -geodtest_SOURCES = geodtest.c +geodtest_SOURCES = geodtest.c geod.c proj_LDADD = libproj.la cs2cs_LDADD = libproj.la @@ -30,7 +30,6 @@ nad2bin_LDADD = libproj.la geod_LDADD = libproj.la multistresstest_LDADD = libproj.la @THREAD_LIB@ test228_LDADD = libproj.la @THREAD_LIB@ -geodtest_LDADD = libproj.la lib_LTLIBRARIES = libproj.la -- cgit v1.2.3 From 195c883b7f76bd5193052f3d5532e73005eadf0c Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Wed, 30 Aug 2017 10:43:48 -0400 Subject: Back out of previous change. Coverage stats are still wrong. geodtest should have included geodesic.c and not geod.c and yet I'm not allowed to use geodesic.c except through the library and then the coverage of geodesic.c seems to be wildly wrong. --- src/Makefile.am | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/Makefile.am b/src/Makefile.am index d215e678..12d11c5f 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -22,7 +22,7 @@ nad2bin_SOURCES = nad2bin.c geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h multistresstest_SOURCES = multistresstest.c test228_SOURCES = test228.c -geodtest_SOURCES = geodtest.c geod.c +geodtest_SOURCES = geodtest.c proj_LDADD = libproj.la cs2cs_LDADD = libproj.la @@ -30,6 +30,7 @@ nad2bin_LDADD = libproj.la geod_LDADD = libproj.la multistresstest_LDADD = libproj.la @THREAD_LIB@ test228_LDADD = libproj.la @THREAD_LIB@ +geodtest_LDADD = libproj.la lib_LTLIBRARIES = libproj.la -- cgit v1.2.3