diff options
| -rw-r--r-- | man/man3/geodesic.3 | 8 | ||||
| -rw-r--r-- | src/apps/geod_set.cpp | 1 | ||||
| -rw-r--r-- | src/geodesic.c | 311 | ||||
| -rw-r--r-- | src/geodesic.h | 141 | ||||
| -rw-r--r-- | src/pipeline.cpp | 2 | ||||
| -rw-r--r-- | src/tests/geodtest.cpp | 213 |
6 files changed, 373 insertions, 303 deletions
diff --git a/man/man3/geodesic.3 b/man/man3/geodesic.3 index c2d3d9df..9cfe9864 100644 --- a/man/man3/geodesic.3 +++ b/man/man3/geodesic.3 @@ -53,9 +53,9 @@ measure angles (latitudes, longitudes, and azimuths) in degrees, unlike the rest of the \fBproj\fR library, which uses radians. The documentation for this library is included in geodesic.h. A formatted version of the documentation is available at -https://geographiclib.sourceforge.io/1.49/C. Detailed documentation of +https://geographiclib.sourceforge.io/1.50/C. Detailed documentation of the interface is given at -https://geographiclib.sourceforge.io/1.49/C/geodesic_8h.html. +https://geographiclib.sourceforge.io/1.50/C/geodesic_8h.html. .SH EXAMPLE The following program reads in lines with the coordinates for two points in decimal degrees (\fIlat1\fR, \fIlon1\fR, \fIlat2\fR, \fIlon2\fR) and @@ -89,9 +89,9 @@ libproj.a \- library of projections and support procedures .SH SEE ALSO Full online documentation for \fBgeodesic(3)\fR, .br -https://geographiclib.sourceforge.io/1.49/C +https://geographiclib.sourceforge.io/1.50/C .br -https://geographiclib.sourceforge.io/1.49/C/geodesic_8h.html +https://geographiclib.sourceforge.io/1.50/C/geodesic_8h.html .PP .B geod(1) .PP diff --git a/src/apps/geod_set.cpp b/src/apps/geod_set.cpp index cd6b5018..ed7edeb9 100644 --- a/src/apps/geod_set.cpp +++ b/src/apps/geod_set.cpp @@ -46,7 +46,6 @@ geod_set(int argc, char **argv) { /* check if line or arc mode */ if (pj_param(nullptr,start, "tlat_1").i) { double del_S; -#undef f phi1 = pj_param(nullptr,start, "rlat_1").f; lam1 = pj_param(nullptr,start, "rlon_1").f; if (pj_param(nullptr,start, "tlat_2").i) { diff --git a/src/geodesic.c b/src/geodesic.c index dcdd4d77..02887e10 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -18,13 +18,15 @@ * * See the comments in geodesic.h for documentation. * - * Copyright (c) Charles Karney (2012-2018) <charles@karney.com> and licensed + * Copyright (c) Charles Karney (2012-2019) <charles@karney.com> and licensed * under the MIT/X11 License. For more information, see * https://geographiclib.sourceforge.io/ */ #include "geodesic.h" #include <math.h> +#include <limits.h> +#include <float.h> #if !defined(HAVE_C99_MATH) #if defined(PROJ_LIB) @@ -65,21 +67,9 @@ static real epsilon, realmin, pi, degree, NaN, static void Init() { if (!init) { -#if defined(__DBL_MANT_DIG__) - digits = __DBL_MANT_DIG__; -#else - digits = 53; -#endif -#if defined(__DBL_EPSILON__) - epsilon = __DBL_EPSILON__; -#else - epsilon = pow(0.5, digits - 1); -#endif -#if defined(__DBL_MIN__) - realmin = __DBL_MIN__; -#else - realmin = pow(0.5, 1022); -#endif + digits = DBL_MANT_DIG; + epsilon = DBL_EPSILON; + realmin = DBL_MIN; #if defined(M_PI) pi = M_PI; #else @@ -98,15 +88,15 @@ static void Init() { tolb = tol0 * tol2; xthresh = 1000 * tol2; degree = pi/180; - #if defined(NAN) - NaN = NAN; - #else +#if defined(NAN) + NaN = NAN; /* NAN is defined in C99 */ +#else { real minus1 = -1; /* cppcheck-suppress wrongmathcall */ NaN = sqrt(minus1); } - #endif +#endif init = 1; } } @@ -122,13 +112,28 @@ enum captype { OUT_ALL = 0x7F80U }; -static real sq(real x) { return x * x; } #if HAVE_C99_MATH +#define hypotx hypot +/* no need to redirect log1px, since it's only used by atanhx */ #define atanhx atanh #define copysignx copysign -#define hypotx hypot #define cbrtx cbrt +#define remainderx remainder +#define remquox remquo #else +/* Replacements for C99 math functions */ + +static real hypotx(real x, real y) { + x = fabs(x); y = fabs(y); + if (x < y) { + x /= y; /* y is nonzero */ + return y * sqrt(1 + x * x); + } else { + y /= (x != 0 ? x : 1); + return x * sqrt(1 + y * y); + } +} + static real log1px(real x) { volatile real y = 1 + x, @@ -143,22 +148,61 @@ static real log1px(real x) { static real atanhx(real x) { real y = fabs(x); /* Enforce odd parity */ y = log1px(2 * y/(1 - y))/2; - return x < 0 ? -y : y; + return x > 0 ? y : (x < 0 ? -y : x); /* atanh(-0.0) = -0.0 */ } static real copysignx(real x, real y) { + /* 1/y trick to get the sign of -0.0 */ return fabs(x) * (y < 0 || (y == 0 && 1/y < 0) ? -1 : 1); } -static real hypotx(real x, real y) -{ return sqrt(x * x + y * y); } - static real cbrtx(real x) { - real y = pow(fabs(x), 1/(real)(3)); /* Return the real cube root */ - return x < 0 ? -y : y; + real y = pow(fabs(x), 1/(real)(3)); /* Return the real cube root */ + return x > 0 ? y : (x < 0 ? -y : x); /* cbrt(-0.0) = -0.0 */ +} + +static real remainderx(real x, real y) { + real z; + y = fabs(y); /* The result doesn't depend on the sign of y */ + z = fmod(x, y); + if (z == 0) + /* This shouldn't be necessary. However, before version 14 (2015), + * Visual Studio had problems dealing with -0.0. Specifically + * VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0 + * python 2.7 on Windows 32-bit machines has the same problem. */ + z = copysignx(z, x); + else if (2 * fabs(z) == y) + z -= fmod(x, 2 * y) - z; /* Implement ties to even */ + else if (2 * fabs(z) > y) + z += (z < 0 ? y : -y); /* Fold remaining cases to (-y/2, y/2) */ + return z; +} + +static real remquox(real x, real y, int* n) { + real z = remainderx(x, y); + if (n) { + real + a = remainderx(x, 2 * y), + b = remainderx(x, 4 * y), + c = remainderx(x, 8 * y); + *n = (a > z ? 1 : (a < z ? -1 : 0)); + *n += (b > a ? 2 : (b < a ? -2 : 0)); + *n += (c > b ? 4 : (c < b ? -4 : 0)); + if (y < 0) *n *= -1; + if (y != 0) { + if (x/y > 0 && *n <= 0) + *n += 8; + else if (x/y < 0 && *n >= 0) + *n -= 8; + } + } + return z; } + #endif +static real sq(real x) { return x * x; } + static real sumx(real u, real v, real* t) { volatile real s = u + v; volatile real up = s - v; @@ -195,23 +239,8 @@ static void norm2(real* sinx, real* cosx) { } static real AngNormalize(real x) { -#if HAVE_C99_MATH - x = remainder(x, (real)(360)); + x = remainderx(x, (real)(360)); return x != -180 ? x : 180; -#else - real y = fmod(x, (real)(360)); -#if defined(_MSC_VER) && _MSC_VER < 1900 - /* - * Before version 14 (2015), Visual Studio had problems dealing - * with -0.0. Specifically - * VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0 - * sincosdx has a similar fix. - * python 2.7 on Windows 32-bit machines has the same problem. - */ - if (x == 0) y = x; -#endif - return y <= -180 ? y + 360 : (y <= 180 ? y : y - 360); -#endif } static real LatFix(real x) @@ -242,16 +271,7 @@ 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)); - /* check for NaN */ - q = r == r ? (int)(floor(r / 90 + (real)(0.5))) : 0; - r -= 90 * q; -#endif + r = remquox(x, (real)(90), &q); /* now abs(r) <= 45 */ r *= degree; /* Possibly could call the gnu extension sincos */ @@ -355,6 +375,11 @@ static void acccopy(const real s[], real t[]); static void accadd(real s[], real y); static real accsum(const real s[], real y); static void accneg(real s[]); +static void accrem(real s[], real y); +static real areareduceA(real area[], real area0, + int crossings, boolx reverse, boolx sign); +static real areareduceB(real area, real area0, + int crossings, boolx reverse, boolx sign); void geod_init(struct geod_geodesic* g, real a, real f) { if (!init) Init(); @@ -698,12 +723,14 @@ real geod_genposition(const struct geod_geodesicline* l, void geod_setdistance(struct geod_geodesicline* l, real s13) { l->s13 = s13; - l->a13 = geod_genposition(l, GEOD_NOFLAGS, l->s13, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); + l->a13 = geod_genposition(l, GEOD_NOFLAGS, l->s13, nullptr, nullptr, nullptr, + nullptr, nullptr, nullptr, nullptr, nullptr); } static void geod_setarc(struct geod_geodesicline* l, real a13) { l->a13 = a13; l->s13 = NaN; - geod_genposition(l, GEOD_ARCMODE, l->a13, nullptr, nullptr, nullptr, &l->s13, nullptr, nullptr, nullptr, nullptr); + geod_genposition(l, GEOD_ARCMODE, l->a13, nullptr, nullptr, nullptr, &l->s13, + nullptr, nullptr, nullptr, nullptr); } void geod_gensetdistance(struct geod_geodesicline* l, @@ -715,7 +742,8 @@ void geod_gensetdistance(struct geod_geodesicline* l, void geod_position(const struct geod_geodesicline* l, real s12, real* plat2, real* plon2, real* pazi2) { - geod_genposition(l, FALSE, s12, plat2, plon2, pazi2, nullptr, nullptr, nullptr, nullptr, nullptr); + geod_genposition(l, FALSE, s12, plat2, plon2, pazi2, + nullptr, nullptr, nullptr, nullptr, nullptr); } real geod_gendirect(const struct geod_geodesic* g, @@ -1134,7 +1162,8 @@ void geod_inverseline(struct geod_geodesicline* l, void geod_inverse(const struct geod_geodesic* g, real lat1, real lon1, real lat2, real lon2, real* ps12, real* pazi1, real* pazi2) { - geod_geninverse(g, lat1, lon1, lat2, lon2, ps12, pazi1, pazi2, nullptr, nullptr, nullptr, nullptr); + geod_geninverse(g, lat1, lon1, lat2, lon2, ps12, pazi1, pazi2, + nullptr, nullptr, nullptr, nullptr); } real SinCosSeries(boolx sinp, real sinx, real cosx, const real c[], int n) { @@ -1297,22 +1326,7 @@ real InverseStart(const struct geod_geodesic* g, boolx shortline = cbet12 >= 0 && sbet12 < (real)(0.5) && cbet2 * lam12 < (real)(0.5); real somg12, comg12, ssig12, csig12; -#if defined(__GNUC__) && __GNUC__ == 4 && \ - (__GNUC_MINOR__ < 6 || defined(__MINGW32__)) - /* Volatile declaration needed to fix inverse cases - * 88.202499451857 0 -88.202499451857 179.981022032992859592 - * 89.262080389218 0 -89.262080389218 179.992207982775375662 - * 89.333123580033 0 -89.333123580032997687 179.99295812360148422 - * which otherwise fail with g++ 4.4.4 x86 -O3 (Linux) - * and g++ 4.4.0 (mingw) and g++ 4.6.1 (tdm mingw). */ - { - volatile real xx1 = sbet2 * cbet1; - volatile real xx2 = cbet2 * sbet1; - sbet12a = xx1 + xx2; - } -#else sbet12a = sbet2 * cbet1 + cbet2 * sbet1; -#endif if (shortline) { real sbetm2 = sq(sbet1 + sbet2), omg12; /* sin((bet1+bet2)/2)^2 @@ -1830,17 +1844,12 @@ 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)); + /* Compute exactly the parity of + int(ceil(lon2 / 360)) - int(ceil(lon1 / 360)) */ + lon1 = remainderx(lon1, (real)(720)); + lon2 = remainderx(lon2, (real)(720)); return ( (lon2 <= 0 && lon2 > -360 ? 1 : 0) - (lon1 <= 0 && lon1 > -360 ? 1 : 0) ); -#else - lon1 = fmod(lon1, (real)(720)); - lon2 = fmod(lon2, (real)(720)); - return ( ((lon2 <= 0 && lon2 > -360) || lon2 > 360 ? 1 : 0) - - ((lon1 <= 0 && lon1 > -360) || lon1 > 360 ? 1 : 0) ); -#endif } void accini(real s[]) { @@ -1876,6 +1885,12 @@ void accneg(real s[]) { s[0] = -s[0]; s[1] = -s[1]; } +void accrem(real s[], real y) { + /* Reduce to [-y/2, y/2]. */ + s[0] = remainderx(s[0], y); + accadd(s, (real)(0)); +} + void geod_polygon_init(struct geod_polygon* p, boolx polylinep) { p->polyline = (polylinep != 0); geod_polygon_clear(p); @@ -1898,7 +1913,8 @@ void geod_polygon_addpoint(const struct geod_geodesic* g, } else { real s12, S12 = 0; /* Initialize S12 to stop Visual Studio warning */ geod_geninverse(g, p->lat, p->lon, lat, lon, - &s12, nullptr, nullptr, nullptr, nullptr, nullptr, p->polyline ? nullptr : &S12); + &s12, nullptr, nullptr, nullptr, nullptr, nullptr, + p->polyline ? nullptr : &S12); accadd(p->P, s12); if (!p->polyline) { accadd(p->A, S12); @@ -1918,7 +1934,8 @@ void geod_polygon_addedge(const struct geod_geodesic* g, real lat = 0, lon = 0, S12 = 0; geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s, &lat, &lon, nullptr, - nullptr, nullptr, nullptr, nullptr, p->polyline ? nullptr : &S12); + nullptr, nullptr, nullptr, nullptr, + p->polyline ? nullptr : &S12); accadd(p->P, s); if (!p->polyline) { accadd(p->A, S12); @@ -1933,8 +1950,7 @@ unsigned geod_polygon_compute(const struct geod_geodesic* g, const struct geod_polygon* p, boolx reverse, boolx sign, real* pA, real* pP) { - real s12, S12, t[2], area0; - int crossings; + real s12, S12, t[2]; if (p->num < 2) { if (pP) *pP = 0; if (!p->polyline && pA) *pA = 0; @@ -1949,27 +1965,9 @@ unsigned geod_polygon_compute(const struct geod_geodesic* g, if (pP) *pP = accsum(p->P, s12); acccopy(p->A, t); accadd(t, S12); - crossings = p->crossings + transit(p->lon, p->lon0); - area0 = 4 * pi * g->c2; - if (crossings & 1) - accadd(t, (t[0] < 0 ? 1 : -1) * area0/2); - /* area is with the clockwise sense. If !reverse convert to - * counter-clockwise convention. */ - if (!reverse) - accneg(t); - /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */ - if (sign) { - if (t[0] > area0/2) - accadd(t, -area0); - else if (t[0] <= -area0/2) - accadd(t, +area0); - } else { - if (t[0] >= area0) - accadd(t, -area0); - else if (t[0] < 0) - accadd(t, +area0); - } - if (pA) *pA = 0 + t[0]; + if (pA) *pA = areareduceA(t, 4 * pi * g->c2, + p->crossings + transit(p->lon, p->lon0), + reverse, sign); return p->num; } @@ -1978,7 +1976,7 @@ unsigned geod_polygon_testpoint(const struct geod_geodesic* g, real lat, real lon, boolx reverse, boolx sign, real* pA, real* pP) { - real perimeter, tempsum, area0; + real perimeter, tempsum; int crossings, i; unsigned num = p->num + 1; if (num == 1) { @@ -1994,7 +1992,8 @@ unsigned geod_polygon_testpoint(const struct geod_geodesic* g, geod_geninverse(g, i == 0 ? p->lat : lat, i == 0 ? p->lon : lon, i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon, - &s12, nullptr, nullptr, nullptr, nullptr, nullptr, p->polyline ? nullptr : &S12); + &s12, nullptr, nullptr, nullptr, nullptr, nullptr, + p->polyline ? nullptr : &S12); perimeter += s12; if (!p->polyline) { tempsum += S12; @@ -2007,26 +2006,7 @@ unsigned geod_polygon_testpoint(const struct geod_geodesic* g, if (p->polyline) return num; - area0 = 4 * pi * g->c2; - if (crossings & 1) - tempsum += (tempsum < 0 ? 1 : -1) * area0/2; - /* area is with the clockwise sense. If !reverse convert to - * counter-clockwise convention. */ - if (!reverse) - tempsum *= -1; - /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */ - if (sign) { - if (tempsum > area0/2) - tempsum -= area0; - else if (tempsum <= -area0/2) - tempsum += area0; - } else { - if (tempsum >= area0) - tempsum -= area0; - else if (tempsum < 0) - tempsum += area0; - } - if (pA) *pA = 0 + tempsum; + if (pA) *pA = areareduceB(tempsum, 4 * pi * g->c2, crossings, reverse, sign); return num; } @@ -2035,7 +2015,7 @@ unsigned geod_polygon_testedge(const struct geod_geodesic* g, real azi, real s, boolx reverse, boolx sign, real* pA, real* pP) { - real perimeter, tempsum, area0; + real perimeter, tempsum; int crossings; unsigned num = p->num + 1; if (num == 1) { /* we don't have a starting point! */ @@ -2053,7 +2033,7 @@ unsigned geod_polygon_testedge(const struct geod_geodesic* g, crossings = p->crossings; { /* Initialization of lat, lon, and S12 is to make CLang static analyzer - happy. */ + * happy. */ real lat = 0, lon = 0, s12, S12 = 0; geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s, &lat, &lon, nullptr, @@ -2067,27 +2047,8 @@ unsigned geod_polygon_testedge(const struct geod_geodesic* g, crossings += transit(lon, p->lon0); } - area0 = 4 * pi * g->c2; - if (crossings & 1) - tempsum += (tempsum < 0 ? 1 : -1) * area0/2; - /* area is with the clockwise sense. If !reverse convert to - * counter-clockwise convention. */ - if (!reverse) - tempsum *= -1; - /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */ - if (sign) { - if (tempsum > area0/2) - tempsum -= area0; - else if (tempsum <= -area0/2) - tempsum += area0; - } else { - if (tempsum >= area0) - tempsum -= area0; - else if (tempsum < 0) - tempsum += area0; - } if (pP) *pP = perimeter; - if (pA) *pA = 0 + tempsum; + if (pA) *pA = areareduceB(tempsum, 4 * pi * g->c2, crossings, reverse, sign); return num; } @@ -2102,4 +2063,52 @@ void geod_polygonarea(const struct geod_geodesic* g, geod_polygon_compute(g, &p, FALSE, TRUE, pA, pP); } +real areareduceA(real area[], real area0, + int crossings, boolx reverse, boolx sign) { + accrem(area, area0); + if (crossings & 1) + accadd(area, (area[0] < 0 ? 1 : -1) * area0/2); + /* area is with the clockwise sense. If !reverse convert to + * counter-clockwise convention. */ + if (!reverse) + accneg(area); + /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */ + if (sign) { + if (area[0] > area0/2) + accadd(area, -area0); + else if (area[0] <= -area0/2) + accadd(area, +area0); + } else { + if (area[0] >= area0) + accadd(area, -area0); + else if (area[0] < 0) + accadd(area, +area0); + } + return 0 + area[0]; +} + +real areareduceB(real area, real area0, + int crossings, boolx reverse, boolx sign) { + area = remainderx(area, area0); + if (crossings & 1) + area += (area < 0 ? 1 : -1) * area0/2; + /* area is with the clockwise sense. If !reverse convert to + * counter-clockwise convention. */ + if (!reverse) + area *= -1; + /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */ + if (sign) { + if (area > area0/2) + area -= area0; + else if (area <= -area0/2) + area += area0; + } else { + if (area >= area0) + area -= area0; + else if (area < 0) + area += area0; + } + return 0 + area; +} + /** @endcond */ diff --git a/src/geodesic.h b/src/geodesic.h index 11484ec7..5d230531 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-2018) <charles@karney.com> and licensed + * Copyright (c) Charles Karney (2012-2019) <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.49. + * <a href="../index.html">GeographicLib</a> 1.50. **********************************************************************/ #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 49 +#define GEODESIC_VERSION_MINOR 50 /** * The patch level of the geodesic library. (This tracks the version of * GeographicLib.) **********************************************************************/ -#define GEODESIC_VERSION_PATCH 3 +#define GEODESIC_VERSION_PATCH 0 /** * Pack the version components into a single integer. Users should not rely on @@ -157,7 +157,7 @@ GEODESIC_VERSION_MINOR, \ GEODESIC_VERSION_PATCH) -#ifndef GEOD_DLL +#if !defined(GEOD_DLL) #if defined(_MSC_VER) #define GEOD_DLL __declspec(dllexport) #elif defined(__GNUC__) @@ -167,7 +167,7 @@ #endif #endif -#ifdef PROJ_RENAME_SYMBOLS +#if defined(PROJ_RENAME_SYMBOLS) #include "proj_symbol_rename.h" #endif @@ -277,8 +277,8 @@ extern "C" { @endcode **********************************************************************/ void GEOD_DLL geod_direct(const struct geod_geodesic* g, - double lat1, double lon1, double azi1, double s12, - double* plat2, double* plon2, double* pazi2); + double lat1, double lon1, double azi1, double s12, + double* plat2, double* plon2, double* pazi2); /** * The general direct geodesic problem. @@ -319,11 +319,12 @@ extern "C" { * what sense the geodesic encircles the ellipsoid. **********************************************************************/ double GEOD_DLL geod_gendirect(const struct geod_geodesic* g, - double lat1, double lon1, double azi1, - unsigned flags, double s12_a12, - double* plat2, double* plon2, double* pazi2, - double* ps12, double* pm12, double* pM12, double* pM21, - double* pS12); + double lat1, double lon1, double azi1, + unsigned flags, double s12_a12, + double* plat2, double* plon2, double* pazi2, + double* ps12, double* pm12, + double* pM12, double* pM21, + double* pS12); /** * Solve the inverse geodesic problem. @@ -364,8 +365,9 @@ extern "C" { @endcode **********************************************************************/ void GEOD_DLL geod_inverse(const struct geod_geodesic* g, - double lat1, double lon1, double lat2, double lon2, - double* ps12, double* pazi1, double* pazi2); + double lat1, double lon1, + double lat2, double lon2, + double* ps12, double* pazi1, double* pazi2); /** * The general inverse geodesic calculation. @@ -395,10 +397,11 @@ extern "C" { * some quantities computed. **********************************************************************/ double GEOD_DLL geod_geninverse(const struct geod_geodesic* g, - double lat1, double lon1, double lat2, double lon2, - double* ps12, double* pazi1, double* pazi2, - double* pm12, double* pM12, double* pM21, - double* pS12); + double lat1, double lon1, + double lat2, double lon2, + double* ps12, double* pazi1, double* pazi2, + double* pm12, double* pM12, double* pM21, + double* pS12); /** * Initialize a geod_geodesicline object. @@ -440,8 +443,9 @@ extern "C" { * NaN). **********************************************************************/ void GEOD_DLL geod_lineinit(struct geod_geodesicline* l, - const struct geod_geodesic* g, - double lat1, double lon1, double azi1, unsigned caps); + const struct geod_geodesic* g, + double lat1, double lon1, double azi1, + unsigned caps); /** * Initialize a geod_geodesicline object in terms of the direct geodesic @@ -465,9 +469,10 @@ extern "C" { * information. **********************************************************************/ void GEOD_DLL geod_directline(struct geod_geodesicline* l, - const struct geod_geodesic* g, - double lat1, double lon1, double azi1, double s12, - unsigned caps); + const struct geod_geodesic* g, + double lat1, double lon1, + double azi1, double s12, + unsigned caps); /** * Initialize a geod_geodesicline object in terms of the direct geodesic @@ -495,10 +500,10 @@ extern "C" { * information. **********************************************************************/ void GEOD_DLL geod_gendirectline(struct geod_geodesicline* l, - const struct geod_geodesic* g, - double lat1, double lon1, double azi1, - unsigned flags, double s12_a12, - unsigned caps); + const struct geod_geodesic* g, + double lat1, double lon1, double azi1, + unsigned flags, double s12_a12, + unsigned caps); /** * Initialize a geod_geodesicline object in terms of the inverse geodesic @@ -521,9 +526,10 @@ extern "C" { * information. **********************************************************************/ void GEOD_DLL geod_inverseline(struct geod_geodesicline* l, - const struct geod_geodesic* g, - double lat1, double lon1, double lat2, double lon2, - unsigned caps); + const struct geod_geodesic* g, + double lat1, double lon1, + double lat2, double lon2, + unsigned caps); /** * Compute the position along a geod_geodesicline. @@ -538,9 +544,9 @@ extern "C" { * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). * * \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 [−180°, 180°]. Any of the - * "return" arguments \e plat2, etc., may be replaced by 0, if you do not + * with \e caps |= GEOD_DISTANCE_IN (or \e caps = 0). The values of \e lon2 + * and \e azi2 returned are in the range [−180°, 180°]. Any of + * the "return" arguments \e plat2, etc., may be replaced by 0, if you do not * need some quantities computed. * * Example, compute way points between JFK and Singapore Changi Airport @@ -571,7 +577,7 @@ extern "C" { @endcode **********************************************************************/ void GEOD_DLL geod_position(const struct geod_geodesicline* l, double s12, - double* plat2, double* plon2, double* pazi2); + double* plat2, double* plon2, double* pazi2); /** * The general position function. @@ -638,11 +644,11 @@ extern "C" { @endcode **********************************************************************/ double GEOD_DLL geod_genposition(const struct geod_geodesicline* l, - unsigned flags, double s12_a12, - double* plat2, double* plon2, double* pazi2, - double* ps12, double* pm12, - double* pM12, double* pM21, - double* pS12); + unsigned flags, double s12_a12, + double* plat2, double* plon2, double* pazi2, + double* ps12, double* pm12, + double* pM12, double* pM21, + double* pS12); /** * Specify position of point 3 in terms of distance. @@ -672,7 +678,7 @@ extern "C" { * been constructed with \e caps |= GEOD_DISTANCE. **********************************************************************/ void GEOD_DLL geod_gensetdistance(struct geod_geodesicline* l, - unsigned flags, double s13_a13); + unsigned flags, double s13_a13); /** * Initialize a geod_polygon object. @@ -721,8 +727,8 @@ extern "C" { * geod_polygon_compute(). **********************************************************************/ void GEOD_DLL geod_polygon_addpoint(const struct geod_geodesic* g, - struct geod_polygon* p, - double lat, double lon); + struct geod_polygon* p, + double lat, double lon); /** * Add an edge to the polygon or polyline. @@ -741,8 +747,8 @@ extern "C" { * new vertex. **********************************************************************/ void GEOD_DLL geod_polygon_addedge(const struct geod_geodesic* g, - struct geod_polygon* p, - double azi, double s); + struct geod_polygon* p, + double azi, double s); /** * Return the results for a polygon. @@ -763,10 +769,12 @@ extern "C" { * * The area and perimeter are accumulated at two times the standard floating * point precision to guard against the loss of accuracy with many-sided - * polygons. Only simple polygons (which are not self-intersecting) are - * allowed. There's no need to "close" the polygon by repeating the first - * vertex. Set \e pA or \e pP to zero, if you do not want the corresponding - * quantity returned. + * polygons. Arbitrarily complex polygons are allowed. In the case of + * self-intersecting polygons the area is accumulated "algebraically", e.g., + * the areas of the 2 loops in a figure-8 polygon will partially cancel. + * There's no need to "close" the polygon by repeating the first vertex. Set + * \e pA or \e pP to zero, if you do not want the corresponding quantity + * returned. * * More points can be added to the polygon after this call. * @@ -788,9 +796,9 @@ extern "C" { @endcode **********************************************************************/ unsigned GEOD_DLL geod_polygon_compute(const struct geod_geodesic* g, - const struct geod_polygon* p, - int reverse, int sign, - double* pA, double* pP); + const struct geod_polygon* p, + int reverse, int sign, + double* pA, double* pP); /** * Return the results assuming a tentative final test point is added; @@ -819,10 +827,10 @@ extern "C" { * \e lat should be in the range [−90°, 90°]. **********************************************************************/ unsigned GEOD_DLL geod_polygon_testpoint(const struct geod_geodesic* g, - const struct geod_polygon* p, - double lat, double lon, - int reverse, int sign, - double* pA, double* pP); + const struct geod_polygon* p, + double lat, double lon, + int reverse, int sign, + double* pA, double* pP); /** * Return the results assuming a tentative final test point is added via an @@ -850,10 +858,10 @@ extern "C" { * @return the number of points. **********************************************************************/ unsigned GEOD_DLL geod_polygon_testedge(const struct geod_geodesic* g, - const struct geod_polygon* p, - double azi, double s, - int reverse, int sign, - double* pA, double* pP); + const struct geod_polygon* p, + double azi, double s, + int reverse, int sign, + double* pA, double* pP); /** * A simple interface for computing the area of a geodesic polygon. @@ -868,10 +876,11 @@ extern "C" { * * \e lats should be in the range [−90°, 90°]. * - * Only simple polygons (which are not self-intersecting) are allowed. - * There's no need to "close" the polygon by repeating the first vertex. The - * area returned is signed with counter-clockwise traversal being treated as - * positive. + * Arbitrarily complex polygons are allowed. In the case self-intersecting + * of polygons the area is accumulated "algebraically", e.g., the areas of + * the 2 loops in a figure-8 polygon will partially cancel. There's no need + * to "close" the polygon by repeating the first vertex. The area returned + * is signed with counter-clockwise traversal being treated as positive. * * Example, compute the area of Antarctica: @code{.c} @@ -888,8 +897,8 @@ extern "C" { @endcode **********************************************************************/ void GEOD_DLL geod_polygonarea(const struct geod_geodesic* g, - double lats[], double lons[], int n, - double* pA, double* pP); + double lats[], double lons[], int n, + double* pA, double* pP); /** * mask values for the \e caps argument to geod_lineinit(). diff --git a/src/pipeline.cpp b/src/pipeline.cpp index afa3b19a..a82fce31 100644 --- a/src/pipeline.cpp +++ b/src/pipeline.cpp @@ -330,7 +330,7 @@ static void set_ellipsoid(PJ *P) { pj_calc_ellipsoid_params (P, P->a, P->es); - geod_init(P->geod, P->a, (1 - sqrt (1 - P->es))); + geod_init(P->geod, P->a, P->es / (1 + sqrt(P->one_es))); /* Re-attach the dangling list */ /* Note: cur will always be non 0 given argv_sentinel presence, */ diff --git a/src/tests/geodtest.cpp b/src/tests/geodtest.cpp index 6b3ea8b2..097682ac 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-2018) <charles@karney.com> and licensed + * Copyright (c) Charles Karney (2015-2019) <charles@karney.com> and licensed * under the MIT/X11 License. For more information, see * https://geographiclib.sourceforge.io/ **********************************************************************/ @@ -20,6 +20,10 @@ # pragma warning (disable: 4706) #endif +#if !defined(__cplusplus) +#define nullptr 0 +#endif + static const double wgs84_a = 6378137, wgs84_f = 1/298.257223563; /* WGS84 */ static int checkEquals(double x, double y, double d) { @@ -30,6 +34,7 @@ static int checkEquals(double x, double y, double d) { } static int checkNaN(double x) { + /* cppcheck-suppress duplicateExpression */ if (x != x) return 0; printf("checkNaN fails: %.7g\n", x); @@ -157,8 +162,7 @@ static int testdirect() { s12 = testcases[i][6]; a12 = testcases[i][7]; m12 = testcases[i][8]; M12 = testcases[i][9]; M21 = testcases[i][10]; S12 = testcases[i][11]; a12a = geod_gendirect(&g, lat1, lon1, azi1, flags, s12, - &lat2a, &lon2a, &azi2a, nullptr, - &m12a, &M12a, &M21a, &S12a); + &lat2a, &lon2a, &azi2a, nullptr, &m12a, &M12a, &M21a, &S12a); result += checkEquals(lat2, lat2a, 1e-13); result += checkEquals(lon2, lon2a, 1e-13); result += checkEquals(azi2, azi2a, 1e-13); @@ -277,13 +281,16 @@ static int GeodSolve6() { int result = 0; geod_init(&g, wgs84_a, wgs84_f); geod_inverse(&g, 88.202499451857, 0, - -88.202499451857, 179.981022032992859592, &s12, nullptr, nullptr); + -88.202499451857, 179.981022032992859592, + &s12, nullptr, nullptr); result += checkEquals(s12, 20003898.214, 0.5e-3); geod_inverse(&g, 89.262080389218, 0, - -89.262080389218, 179.992207982775375662, &s12, nullptr, nullptr); + -89.262080389218, 179.992207982775375662, + &s12, nullptr, nullptr); result += checkEquals(s12, 20003925.854, 0.5e-3); geod_inverse(&g, 89.333123580033, 0, - -89.333123580032997687, 179.99295812360148422, &s12, nullptr, nullptr); + -89.333123580032997687, 179.99295812360148422, + &s12, nullptr, nullptr); result += checkEquals(s12, 20003926.881, 0.5e-3); return result; } @@ -295,7 +302,8 @@ static int GeodSolve9() { int result = 0; geod_init(&g, wgs84_a, wgs84_f); geod_inverse(&g, 56.320923501171, 0, - -56.320923501171, 179.664747671772880215, &s12, nullptr, nullptr); + -56.320923501171, 179.664747671772880215, + &s12, nullptr, nullptr); result += checkEquals(s12, 19993558.287, 0.5e-3); return result; } @@ -308,7 +316,8 @@ static int GeodSolve10() { int result = 0; geod_init(&g, wgs84_a, wgs84_f); geod_inverse(&g, 52.784459512564, 0, - -52.784459512563990912, 179.634407464943777557, &s12, nullptr, nullptr); + -52.784459512563990912, 179.634407464943777557, + &s12, nullptr, nullptr); result += checkEquals(s12, 19991596.095, 0.5e-3); return result; } @@ -321,7 +330,8 @@ static int GeodSolve11() { int result = 0; geod_init(&g, wgs84_a, wgs84_f); geod_inverse(&g, 48.522876735459, 0, - -48.52287673545898293, 179.599720456223079643, &s12, nullptr, nullptr); + -48.52287673545898293, 179.599720456223079643, + &s12, nullptr, nullptr); result += checkEquals(s12, 19989144.774, 0.5e-3); return result; } @@ -366,8 +376,8 @@ static int GeodSolve15() { struct geod_geodesic g; int result = 0; geod_init(&g, 6.4e6, -1/150.0); - geod_gendirect(&g, 1, 2, 3, 0, 4, - nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, &S12); + geod_gendirect(&g, 1, 2, 3, 0, 4, nullptr, nullptr, nullptr, + nullptr, nullptr, nullptr, nullptr, &S12); result += checkEquals(S12, 23700, 0.5); return result; } @@ -380,13 +390,14 @@ static int GeodSolve17() { int result = 0; unsigned flags = GEOD_LONG_UNROLL; geod_init(&g, wgs84_a, wgs84_f); - geod_gendirect(&g, 40, -75, -10, flags, 2e7, - &lat2, &lon2, &azi2, nullptr, nullptr, nullptr, nullptr, nullptr); + geod_gendirect(&g, 40, -75, -10, flags, 2e7, &lat2, &lon2, &azi2, + nullptr, nullptr, nullptr, nullptr, nullptr); result += checkEquals(lat2, -39, 1); result += checkEquals(lon2, -254, 1); result += checkEquals(azi2, -170, 1); geod_lineinit(&l, &g, 40, -75, -10, 0); - geod_genposition(&l, flags, 2e7, &lat2, &lon2, &azi2, nullptr, nullptr, nullptr, nullptr, nullptr); + geod_genposition(&l, flags, 2e7, &lat2, &lon2, &azi2, + nullptr, nullptr, nullptr, nullptr, nullptr); result += checkEquals(lat2, -39, 1); result += checkEquals(lon2, -254, 1); result += checkEquals(azi2, -170, 1); @@ -407,7 +418,8 @@ static int GeodSolve26() { struct geod_geodesic g; int result = 0; geod_init(&g, 6.4e6, 0); - geod_geninverse(&g, 1, 2, 3, 4, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, &S12); + geod_geninverse(&g, 1, 2, 3, 4, nullptr, nullptr, nullptr, + nullptr, nullptr, nullptr, &S12); result += checkEquals(S12, 49911046115.0, 0.5); return result; } @@ -419,7 +431,8 @@ static int GeodSolve28() { struct geod_geodesic g; int result = 0; geod_init(&g, 6.4e6, 0.1); - a12 = geod_gendirect(&g, 1, 2, 10, 0, 5e6, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); + a12 = geod_gendirect(&g, 1, 2, 10, 0, 5e6, nullptr, nullptr, nullptr, + nullptr, nullptr, nullptr, nullptr, nullptr); result += checkEquals(a12, 48.55570690, 0.5e-8); return result; } @@ -527,12 +540,14 @@ static int GeodSolve61() { unsigned flags = GEOD_LONG_UNROLL; geod_init(&g, wgs84_a, wgs84_f); geod_gendirect(&g, 45, 0, -0.000000000000000003, flags, 1e7, - &lat2, &lon2, &azi2, nullptr, nullptr, nullptr, nullptr, nullptr); + &lat2, &lon2, &azi2, + nullptr, nullptr, nullptr, nullptr, nullptr); result += checkEquals(lat2, 45.30632, 0.5e-5); result += checkEquals(lon2, -180, 0.5e-5); result += checkEquals(fabs(azi2), 180, 0.5e-5); geod_inverseline(&l, &g, 45, 0, 80, -0.000000000000000003, 0); - geod_genposition(&l, flags, 1e7, &lat2, &lon2, &azi2, nullptr, nullptr, nullptr, nullptr, nullptr); + geod_genposition(&l, flags, 1e7, &lat2, &lon2, &azi2, + nullptr, nullptr, nullptr, nullptr, nullptr); result += checkEquals(lat2, 45.30632, 0.5e-5); result += checkEquals(lon2, -180, 0.5e-5); result += checkEquals(fabs(azi2), 180, 0.5e-5); @@ -577,7 +592,7 @@ static int GeodSolve65() { static int GeodSolve67() { /* Check for InverseLine if line is slightly west of S and that s13 is - correctly set. */ + * correctly set. */ double lat2, lon2, azi2; struct geod_geodesic g; struct geod_geodesicline l; @@ -585,11 +600,13 @@ static int GeodSolve67() { unsigned flags = GEOD_LONG_UNROLL; geod_init(&g, wgs84_a, wgs84_f); geod_inverseline(&l, &g, -5, -0.000000000000002, -10, 180, 0); - geod_genposition(&l, flags, 2e7, &lat2, &lon2, &azi2, nullptr, nullptr, nullptr, nullptr, nullptr); + geod_genposition(&l, flags, 2e7, &lat2, &lon2, &azi2, + nullptr, nullptr, nullptr, nullptr, nullptr); result += checkEquals(lat2, 4.96445, 0.5e-5); result += checkEquals(lon2, -180.00000, 0.5e-5); result += checkEquals(azi2, -0.00000, 0.5e-5); - geod_genposition(&l, flags, 0.5 * l.s13, &lat2, &lon2, &azi2, nullptr, nullptr, nullptr, nullptr, nullptr); + geod_genposition(&l, flags, 0.5 * l.s13, &lat2, &lon2, &azi2, + nullptr, nullptr, nullptr, nullptr, nullptr); result += checkEquals(lat2, -87.52461, 0.5e-5); result += checkEquals(lon2, -0.00000, 0.5e-5); result += checkEquals(azi2, -180.00000, 0.5e-5); @@ -614,7 +631,9 @@ static int GeodSolve71() { static int GeodSolve73() { /* Check for backwards from the pole bug reported by Anon on 2016-02-13. * This only affected the Java implementation. It was introduced in Java - * version 1.44 and fixed in 1.46-SNAPSHOT on 2016-01-17. */ + * version 1.44 and fixed in 1.46-SNAPSHOT on 2016-01-17. + * Also the + sign on azi2 is a check on the normalizing of azimuths + * (converting -0.0 to +0.0). */ double lat2, lon2, azi2; struct geod_geodesic g; int result = 0; @@ -652,7 +671,7 @@ static void polylength(const struct geod_geodesic* g, static int GeodSolve74() { /* Check fix for inaccurate areas, bug introduced in v1.46, fixed - 2015-10-16. */ + * 2015-10-16. */ double a12, s12, azi1, azi2, m12, M12, M21, S12; struct geod_geodesic g; int result = 0; @@ -672,7 +691,7 @@ static int GeodSolve74() { static int GeodSolve76() { /* The distance from Wellington and Salamanca (a classic failure of - Vincenty) */ + * Vincenty) */ double azi1, azi2, s12; struct geod_geodesic g; int result = 0; @@ -700,19 +719,24 @@ static int GeodSolve78() { static int GeodSolve80() { /* Some tests to add code coverage: computing scale in special cases + zero - length geodesic (includes GeodSolve80 - GeodSolve83) + using an incapable - line. */ + * length geodesic (includes GeodSolve80 - GeodSolve83) + using an incapable + * line. */ double a12, s12, azi1, azi2, m12, M12, M21, S12; struct geod_geodesic g; struct geod_geodesicline l; int result = 0; geod_init(&g, wgs84_a, wgs84_f); - geod_geninverse(&g, 0, 0, 0, 90, nullptr, nullptr, nullptr, nullptr, &M12, &M21, nullptr); + + geod_geninverse(&g, 0, 0, 0, 90, nullptr, nullptr, nullptr, + nullptr, &M12, &M21, nullptr); result += checkEquals(M12, -0.00528427534, 0.5e-10); result += checkEquals(M21, -0.00528427534, 0.5e-10); - geod_geninverse(&g, 0, 0, 1e-6, 1e-6, nullptr, nullptr, nullptr, nullptr, &M12, &M21, nullptr); + + geod_geninverse(&g, 0, 0, 1e-6, 1e-6, nullptr, nullptr, nullptr, + nullptr, &M12, &M21, nullptr); result += checkEquals(M12, 1, 0.5e-10); result += checkEquals(M21, 1, 0.5e-10); + a12 = geod_geninverse(&g, 20.001, 0, 20.001, 0, &s12, &azi1, &azi2, &m12, &M12, &M21, &S12); result += checkEquals(a12, 0, 1e-13); @@ -723,6 +747,7 @@ static int GeodSolve80() { result += checkEquals(M12, 1, 1e-15); result += checkEquals(M21, 1, 1e-15); result += checkEquals(S12, 0, 1e-10); + a12 = geod_geninverse(&g, 90, 0, 90, 180, &s12, &azi1, &azi2, &m12, &M12, &M21, &S12); result += checkEquals(a12, 0, 1e-13); @@ -732,14 +757,70 @@ static int GeodSolve80() { result += checkEquals(m12, 0, 1e-8); result += checkEquals(M12, 1, 1e-15); result += checkEquals(M21, 1, 1e-15); - result += checkEquals(S12, 127516405431022, 0.5); + result += checkEquals(S12, 127516405431022.0, 0.5); + /* An incapable line which can't take distance as input */ geod_lineinit(&l, &g, 1, 2, 90, GEOD_LATITUDE); - a12 = geod_genposition(&l, 0, 1000, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); + a12 = geod_genposition(&l, 0, 1000, nullptr, nullptr, nullptr, + nullptr, nullptr, nullptr, nullptr, nullptr); result += checkNaN(a12); return result; } +static int GeodSolve84() { + /* Tests for python implementation to check fix for range errors with + * {fmod,sin,cos}(inf) (includes GeodSolve84 - GeodSolve86). */ + + double lat2, lon2, azi2, inf, nan; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + { + /* a round about way to set inf = 0 */ + geod_direct(&g, 0, 0, 90, 0, &inf, nullptr, nullptr); + /* so that this doesn't give a compiler time error on Windows */ + inf = 1.0/inf; + } + { + double minus1 = -1; + /* cppcheck-suppress wrongmathcall */ + nan = sqrt(minus1); + } + geod_direct(&g, 0, 0, 90, inf, &lat2, &lon2, &azi2); + result += checkNaN(lat2); + result += checkNaN(lon2); + result += checkNaN(azi2); + geod_direct(&g, 0, 0, 90, nan, &lat2, &lon2, &azi2); + result += checkNaN(lat2); + result += checkNaN(lon2); + result += checkNaN(azi2); + geod_direct(&g, 0, 0, inf, 1000, &lat2, &lon2, &azi2); + result += checkNaN(lat2); + result += checkNaN(lon2); + result += checkNaN(azi2); + geod_direct(&g, 0, 0, nan, 1000, &lat2, &lon2, &azi2); + result += checkNaN(lat2); + result += checkNaN(lon2); + result += checkNaN(azi2); + geod_direct(&g, 0, inf, 90, 1000, &lat2, &lon2, &azi2); + result += lat2 == 0 ? 0 : 1; + result += checkNaN(lon2); + result += azi2 == 90 ? 0 : 1; + geod_direct(&g, 0, nan, 90, 1000, &lat2, &lon2, &azi2); + result += lat2 == 0 ? 0 : 1; + result += checkNaN(lon2); + result += azi2 == 90 ? 0 : 1; + geod_direct(&g, inf, 0, 90, 1000, &lat2, &lon2, &azi2); + result += checkNaN(lat2); + result += checkNaN(lon2); + result += checkNaN(azi2); + geod_direct(&g, nan, 0, 90, 1000, &lat2, &lon2, &azi2); + result += checkNaN(lat2); + result += checkNaN(lon2); + result += checkNaN(azi2); + 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}}; @@ -841,7 +922,7 @@ static int Planimeter13() { static int Planimeter15() { /* Coverage tests, includes Planimeter15 - Planimeter18 (combinations of - reverse and sign) + calls to testpoint, testedge, geod_polygonarea. */ + * reverse and sign) + calls to testpoint, testedge, geod_polygonarea. */ struct geod_geodesic g; struct geod_polygon p; double lat[] = {2, 1, 3}, lon[] = {1, 2, 3}; @@ -886,7 +967,7 @@ static int Planimeter15() { static int Planimeter19() { /* Coverage tests, includes Planimeter19 - Planimeter20 (degenerate - polygons) + extra cases. */ + * polygons) + extra cases. */ struct geod_geodesic g; struct geod_polygon p; double area, perim; @@ -916,16 +997,17 @@ static int Planimeter19() { geod_polygon_addpoint(&g, &p, 1, 1); geod_polygon_compute(&g, &p, 0, 1, nullptr, &perim); result += perim == 0 ? 0 : 1; + geod_polygon_addpoint(&g, &p, 1, 1); + geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, nullptr, &perim); + result += checkEquals(perim, 1000, 1e-10); + geod_polygon_testpoint(&g, &p, 2, 2, 0, 1, nullptr, &perim); + result += checkEquals(perim, 156876.149, 0.5e-3); return result; } static int Planimeter21() { /* Some test to add code coverage: multiple circlings of pole (includes - Planimeter21 - Planimeter28) + invocations via testpoint and testedge. - Some of the results for i = 4 in the loop are wrong because we don't - reduce the area to the allowed range correctly. However these cases are - not "simple" polygons, so we'll defer fixing the problem for now. - */ + * Planimeter21 - Planimeter28) + invocations via testpoint and testedge. */ struct geod_geodesic g; struct geod_polygon p; double area, lat = 45, @@ -945,35 +1027,35 @@ static int Planimeter21() { geod_polygon_addpoint(&g, &p, lat, 60); geod_polygon_addpoint(&g, &p, lat, 180); geod_polygon_testpoint(&g, &p, lat, -60, 0, 1, &area, nullptr); - if (i != 4) result += checkEquals(area, i*r, 0.5); + result += checkEquals(area, i*r, 0.5); geod_polygon_testpoint(&g, &p, lat, -60, 0, 0, &area, nullptr); - if (i != 4) result += checkEquals(area, i*r, 0.5); + result += checkEquals(area, i*r, 0.5); geod_polygon_testpoint(&g, &p, lat, -60, 1, 1, &area, nullptr); - if (i != 4) result += checkEquals(area, -i*r, 0.5); + result += checkEquals(area, -i*r, 0.5); geod_polygon_testpoint(&g, &p, lat, -60, 1, 0, &area, nullptr); result += checkEquals(area, -i*r + a0, 0.5); geod_polygon_testedge(&g, &p, a, s, 0, 1, &area, nullptr); - if (i != 4) result += checkEquals(area, i*r, 0.5); + result += checkEquals(area, i*r, 0.5); geod_polygon_testedge(&g, &p, a, s, 0, 0, &area, nullptr); - if (i != 4) result += checkEquals(area, i*r, 0.5); + result += checkEquals(area, i*r, 0.5); geod_polygon_testedge(&g, &p, a, s, 1, 1, &area, nullptr); - if (i != 4) result += checkEquals(area, -i*r, 0.5); + result += checkEquals(area, -i*r, 0.5); geod_polygon_testedge(&g, &p, a, s, 1, 0, &area, nullptr); result += checkEquals(area, -i*r + a0, 0.5); geod_polygon_addpoint(&g, &p, lat, -60); geod_polygon_compute(&g, &p, 0, 1, &area, nullptr); - if (i != 4) result += checkEquals(area, i*r, 0.5); + result += checkEquals(area, i*r, 0.5); geod_polygon_compute(&g, &p, 0, 0, &area, nullptr); - if (i != 4) result += checkEquals(area, i*r, 0.5); + result += checkEquals(area, i*r, 0.5); geod_polygon_compute(&g, &p, 1, 1, &area, nullptr); - if (i != 4) result += checkEquals(area, -i*r, 0.5); + result += checkEquals(area, -i*r, 0.5); geod_polygon_compute(&g, &p, 1, 0, &area, nullptr); result += checkEquals(area, -i*r + a0, 0.5); } return result; } -static int AddEdge1() { +static int Planimeter29() { /* Check fix to transitdirect vs transit zero handling inconsistency */ struct geod_geodesic g; struct geod_polygon p; @@ -986,41 +1068,12 @@ static int AddEdge1() { geod_polygon_addedge(&g, &p, 0, 1000); geod_polygon_addedge(&g, &p, -90, 1000); geod_polygon_compute(&g, &p, 0, 1, &area, nullptr); + /* The area should be 1e6. Prior to the fix it was 1e6 - A/2, where + * A = ellipsoid area. */ result += checkEquals(area, 1000000.0, 0.01); return result; } -static int EmptyPoly() { - struct geod_geodesic g; - struct geod_polygon p; - double perim, area; - int result = 0; - geod_init(&g, wgs84_a, wgs84_f); - geod_polygon_init(&p, 0); - geod_polygon_testpoint(&g, &p, 1, 1, 0, 1, &area, &perim); - result += area == 0 ? 0 : 1; - result += perim == 0 ? 0 : 1; - geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, &area, &perim); - result += checkNaN(area); - result += checkNaN(perim); - geod_polygon_compute(&g, &p, 0, 1, &area, &perim); - result += area == 0 ? 0 : 1; - result += perim == 0 ? 0 : 1; - geod_polygon_init(&p, 1); - geod_polygon_testpoint(&g, &p, 1, 1, 0, 1, nullptr, &perim); - result += perim == 0 ? 0 : 1; - geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, nullptr, &perim); - result += checkNaN(perim); - geod_polygon_compute(&g, &p, 0, 1, nullptr, &perim); - result += perim == 0 ? 0 : 1; - geod_polygon_addpoint(&g, &p, 1, 1); - geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, nullptr, &perim); - result += checkEquals(perim, 1000, 1e-10); - geod_polygon_testpoint(&g, &p, 2, 2, 0, 1, nullptr, &perim); - result += checkEquals(perim, 156876.149, 0.5e-3); - return result; -} - int main() { int n = 0, i; if ((i = testinverse())) {++n; printf("testinverse fail: %d\n", i);} @@ -1053,6 +1106,7 @@ int main() { if ((i = GeodSolve76())) {++n; printf("GeodSolve76 fail: %d\n", i);} 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 = 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);} @@ -1061,8 +1115,7 @@ int main() { if ((i = Planimeter15())) {++n; printf("Planimeter15 fail: %d\n", i);} if ((i = Planimeter19())) {++n; printf("Planimeter19 fail: %d\n", i);} if ((i = Planimeter21())) {++n; printf("Planimeter21 fail: %d\n", i);} - if ((i = AddEdge1())) {++n; printf("AddEdge1 fail: %d\n", i);} - if ((i = EmptyPoly())) {++n; printf("EmptyPoly fail: %d\n", i);} + if ((i = Planimeter29())) {++n; printf("Planimeter29 fail: %d\n", i);} return n; } |