diff options
| -rw-r--r-- | CMakeLists.txt | 6 | ||||
| -rw-r--r-- | configure.ac | 2 | ||||
| -rw-r--r-- | man/man3/geodesic.3 | 12 | ||||
| -rw-r--r-- | src/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | src/Makefile.am | 2 | ||||
| -rw-r--r-- | src/bin_geodtest.cmake | 12 | ||||
| -rw-r--r-- | src/geodesic.c | 275 | ||||
| -rw-r--r-- | src/geodesic.h | 479 | ||||
| -rw-r--r-- | src/geodtest.c | 768 | ||||
| -rw-r--r-- | src/pj_release.c | 2 | ||||
| -rw-r--r-- | src/proj_api.h | 2 |
11 files changed, 1282 insertions, 279 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index a6c281a3..9ad9417e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,9 +31,9 @@ colormsg(_HIBLUE_ "Configuring PROJ:") #PROJ version information ################################################################################# include(Proj4Version) -proj_version(MAJOR 4 MINOR 9 PATCH 2) -set(PROJ_API_VERSION "10") -set(PROJ_BUILD_VERSION "10.0.0") +proj_version(MAJOR 4 MINOR 9 PATCH 3) +set(PROJ_API_VERSION "11") +set(PROJ_BUILD_VERSION "11.0.0") ################################################################################# # Build features and variants diff --git a/configure.ac b/configure.ac index 16343ba1..7c77490e 100644 --- a/configure.ac +++ b/configure.ac @@ -1,7 +1,7 @@ dnl Process this file with autoconf to produce a configure script. AC_PREREQ(2.59) -AC_INIT([PROJ.4 Projections], 4.9.2, [warmerdam@pobox.com], proj) +AC_INIT([PROJ.4 Projections], 4.9.3, [warmerdam@pobox.com], proj) AC_CONFIG_MACRO_DIR([m4]) AC_LANG(C) diff --git a/man/man3/geodesic.3 b/man/man3/geodesic.3 index bf0d764d..450f75f2 100644 --- a/man/man3/geodesic.3 +++ b/man/man3/geodesic.3 @@ -38,7 +38,7 @@ 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 -http://geographiclib.sf.net/1.44/C +http://geographiclib.sf.net/1.46/C .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 @@ -57,7 +57,7 @@ int main() { struct geod_geodesic g; geod_init(&g, a, f); - while (scanf("%lf %lf %lf %lf\en", + while (scanf("%lf %lf %lf %lf", &lat1, &lon1, &lat2, &lon2) == 4) { geod_inverse(&g, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2); @@ -72,7 +72,7 @@ libproj.a \- library of projections and support procedures .SH SEE ALSO Full online documentation for \fBgeodesic(3)\fR, .br -http://geographiclib.sf.net/1.44/C +http://geographiclib.sf.net/1.46/C .PP .B geod(1) .PP @@ -90,8 +90,12 @@ DOI: http://dx.doi.org/10.1007/s00190-012-0578-z .br http://geographiclib.sf.net/geod-addenda.html .PP -The \fIonline geodesic bibliography\fR, +\fIA geodesic bibliography\fR, .br http://geographiclib.sf.net/geodesic-papers/biblio.html +.PP +The Wikipedia page, \fIGeodesics on an ellipsoid\fR, +.br +https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid .SH HOME PAGE https://github.com/OSGeo/proj.4/wiki diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index aba5b4c4..8d7e7d1c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -32,6 +32,7 @@ endif(BUILD_PROJ) if(BUILD_GEOD) include(bin_geod.cmake) + include(bin_geodtest.cmake) endif(BUILD_GEOD) if(BUILD_NAD2BIN) diff --git a/src/Makefile.am b/src/Makefile.am index 83c5f34e..73e3bdb1 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -29,7 +29,7 @@ test228_LDADD = libproj.la @THREAD_LIB@ lib_LTLIBRARIES = libproj.la -libproj_la_LDFLAGS = -no-undefined -version-info 10:0:0 +libproj_la_LDFLAGS = -no-undefined -version-info 11:0:0 libproj_la_SOURCES = \ pj_list.h \ diff --git a/src/bin_geodtest.cmake b/src/bin_geodtest.cmake new file mode 100644 index 00000000..1604385c --- /dev/null +++ b/src/bin_geodtest.cmake @@ -0,0 +1,12 @@ +set(GEODTEST_SRC geodtest.c ) +set(GEODTEST_INCLUDE geod_interface.h) + +source_group("Source Files\\Bin" FILES ${GEODTEST_SRC} ${GEODTEST_INCLUDE}) + +#Executable +add_executable(geodtest ${GEODTEST_SRC} ${GEODTEST_INCLUDE}) +target_link_libraries(geodtest ${PROJ_LIBRARIES}) +# Do not install + +# Instead run as a test +add_test (NAME geodesic-test COMMAND geodtest) diff --git a/src/geodesic.c b/src/geodesic.c index 2dee45e0..8d9c9285 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -14,11 +14,11 @@ * Algorithms for geodesics, * J. Geodesy <b>87</b>, 43--55 (2013); * https://dx.doi.org/10.1007/s00190-012-0578-z - * Addenda: http://geographiclib.sf.net/geod-addenda.html + * Addenda: http://geographiclib.sourceforge.net/geod-addenda.html * * See the comments in geodesic.h for documentation. * - * Copyright (c) Charles Karney (2012-2015) <charles@karney.com> and licensed + * Copyright (c) Charles Karney (2012-2016) <charles@karney.com> and licensed * under the MIT/X11 License. For more information, see * http://geographiclib.sourceforge.net/ */ @@ -119,6 +119,10 @@ static real atanhx(real x) { return x < 0 ? -y : y; } +static real copysignx(real x, real y) { + 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); } @@ -133,7 +137,7 @@ static real sumx(real u, real v, real* t) { volatile real vpp = s - up; up -= u; vpp -= v; - *t = -(up + vpp); + if (t) *t = -(up + vpp); /* error-free sum: * u + v = s + t * = round(u + v) + t */ @@ -146,11 +150,12 @@ static real polyval(int N, const real p[], real x) { return y; } -static real minx(real x, real y) -{ return x < y ? x : y; } +/* mimic C++ std::min and std::max */ +static real minx(real a, real b) +{ return (b < a) ? b : a; } -static real maxx(real x, real y) -{ return x > y ? x : y; } +static real maxx(real a, real b) +{ return (a < b) ? b : a; } static void swapx(real* x, real* y) { real t = *x; *x = *y; *y = t; } @@ -169,7 +174,7 @@ static real AngNormalize(real x) { static real LatFix(real x) { return fabs(x) > 90 ? NaN : x; } -static real AngDiff(real x, real y) { +static real AngDiff(real x, real y, real* e) { real t, d = - AngNormalize(sumx(AngNormalize(x), AngNormalize(-y), &t)); /* Here y - x = d - t (mod 360), exactly, where d is in (-180,180] and * abs(t) <= eps (eps = 2^-45 for doubles). The only case where the @@ -177,15 +182,16 @@ static real AngDiff(real x, real y) { * and t < 0. The case, d = -180 + eps, t = eps, can't happen, since * sum would have returned the exact result in such a case (i.e., given t * = 0). */ - return (d == 180 && t < 0 ? -180 : d) - t; + return sumx(d == 180 && t < 0 ? -180 : d, -t, e); } static real AngRound(real x) { const real z = 1/(real)(16); + if (x == 0) return 0; volatile real y = fabs(x); /* The compiler mustn't "simplify" z - (z - y) to y */ y = y < z ? z - (z - y) : y; - return x < 0 ? 0 - y : y; + return x < 0 ? -y : y; } static void sincosdx(real x, real* sinx, real* cosx) { @@ -249,7 +255,7 @@ static real Astroid(real x, real y); static real InverseStart(const struct geod_geodesic* g, real sbet1, real cbet1, real dn1, real sbet2, real cbet2, real dn2, - real lam12, + real lam12, real slam12, real clam12, real* psalp1, real* pcalp1, /* Only updated if return val >= 0 */ real* psalp2, real* pcalp2, @@ -261,11 +267,13 @@ static real Lambda12(const struct geod_geodesic* g, real sbet1, real cbet1, real dn1, real sbet2, real cbet2, real dn2, real salp1, real calp1, + real slam120, real clam120, real* psalp2, real* pcalp2, real* psig12, real* pssig1, real* pcsig1, real* pssig2, real* pcsig2, - real* peps, real* pdomg12, + real* peps, + real* psomg12, real* pcomg12, boolx diffp, real* pdlam12, /* Scratch area of the right size */ real Ca[]); @@ -288,7 +296,7 @@ static void accneg(real s[]); void geod_init(struct geod_geodesic* g, real a, real f) { if (!init) Init(); g->a = a; - g->f = f <= 1 ? f : 1/f; + g->f = f; g->f1 = 1 - g->f; g->e2 = g->f * (2 - g->f); g->ep2 = g->e2 / sq(g->f1); /* e2 / (1 - e2) */ @@ -315,9 +323,11 @@ void geod_init(struct geod_geodesic* g, real a, real f) { C4coeff(g); } -void geod_lineinit(struct geod_geodesicline* l, - const struct geod_geodesic* g, - real lat1, real lon1, real azi1, unsigned caps) { +static void geod_lineinit_int(struct geod_geodesicline* l, + const struct geod_geodesic* g, + real lat1, real lon1, + real azi1, real salp1, real calp1, + unsigned caps) { real cbet1, sbet1, eps; l->a = g->a; l->f = g->f; @@ -331,9 +341,9 @@ void geod_lineinit(struct geod_geodesicline* l, l->lat1 = LatFix(lat1); l->lon1 = lon1; - l->azi1 = AngNormalize(azi1); - /* Guard against underflow in salp0 */ - sincosdx(AngRound(l->azi1), &l->salp1, &l->calp1); + l->azi1 = azi1; + l->salp1 = salp1; + l->calp1 = calp1; sincosdx(AngRound(l->lat1), &sbet1, &cbet1); sbet1 *= l->f1; /* Ensure cbet1 = +epsilon at poles */ @@ -396,6 +406,34 @@ void geod_lineinit(struct geod_geodesicline* l, l->A4 = sq(l->a) * l->calp0 * l->salp0 * g->e2; l->B41 = SinCosSeries(FALSE, l->ssig1, l->csig1, l->C4a, nC4); } + + l->a13 = l->s13 = NaN; +} + +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; + /* Guard against underflow in salp0 */ + sincosdx(AngRound(azi1), &salp1, &calp1); + geod_lineinit_int(l, g, lat1, lon1, azi1, salp1, calp1, caps); +} + +void geod_gendirectline(struct geod_geodesicline* l, + const struct geod_geodesic* g, + real lat1, real lon1, real azi1, + unsigned flags, real a12_s12, + unsigned caps) { + geod_lineinit(l, g, lat1, lon1, azi1, caps); + geod_gensetdistance(l, flags, a12_s12); +} + +void geod_directline(struct geod_geodesicline* l, + const struct geod_geodesic* g, + real lat1, real lon1, real azi1, + real s12, unsigned caps) { + geod_gendirectline(l, g, lat1, lon1, azi1, GEOD_NOFLAGS, s12, caps); } real geod_genposition(const struct geod_geodesicline* l, @@ -421,7 +459,7 @@ real geod_genposition(const struct geod_geodesicline* l, outmask &= l->caps & OUT_ALL; if (!( TRUE /*Init()*/ && - (flags & GEOD_ARCMODE || (l->caps & GEOD_DISTANCE_IN & OUT_ALL)) )) + (flags & GEOD_ARCMODE || (l->caps & (GEOD_DISTANCE_IN & OUT_ALL))) )) /* Uninitialized or impossible distance calculation requested */ return NaN; @@ -464,10 +502,9 @@ real geod_genposition(const struct geod_geodesicline* l, * 1/20 5376 146e3 10e3 * 1/10 829e3 22e6 1.5e6 * 1/5 157e6 3.8e9 280e6 */ - real - ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12, - csig2 = l->csig1 * csig12 - l->ssig1 * ssig12, - serr; + real serr; + ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12; + csig2 = l->csig1 * csig12 - l->ssig1 * ssig12; B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1); serr = (1 + l->A1m1) * (sig12 + (B12 - l->B11)) - s12_a12 / l->b; sig12 = sig12 - serr / sqrt(1 + l->k2 * sq(ssig2)); @@ -499,7 +536,7 @@ real geod_genposition(const struct geod_geodesicline* l, s12 = flags & GEOD_ARCMODE ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12; if (outmask & GEOD_LONGITUDE) { - int E = l->salp0 < 0 ? -1 : 1; /* east or west going? */ + real E = copysignx(1, l->salp0); /* east or west going? */ /* tan(omg2) = sin(alp0) * tan(sig2) */ somg2 = l->salp0 * ssig2; comg2 = csig2; /* No need to normalize */ /* omg12 = omg2 - omg1 */ @@ -548,14 +585,6 @@ real geod_genposition(const struct geod_geodesicline* l, /* alp12 = alp2 - alp1, used in atan2 so no need to normalize */ salp12 = salp2 * l->calp1 - calp2 * l->salp1; calp12 = calp2 * l->calp1 + salp2 * l->salp1; - /* The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz - * salp12 = -0 and alp12 = -180. However this depends on the sign being - * attached to 0 correctly. The following ensures the correct - * behavior. */ - if (salp12 == 0 && calp12 < 0) { - salp12 = tiny * l->calp1; - calp12 = -1; - } } else { /* tan(alp) = tan(alp0) * sec(sig) * tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1) @@ -593,6 +622,21 @@ real geod_genposition(const struct geod_geodesicline* l, return flags & GEOD_ARCMODE ? s12_a12 : sig12 / degree; } +void geod_setdistance(struct geod_geodesicline* l, real s13) { + l->s13 = s13; + l->a13 = geod_genposition(l, GEOD_NOFLAGS, l->s13, 0, 0, 0, 0, 0, 0, 0, 0); +} + +static void geod_setarc(struct geod_geodesicline* l, real a13) { + l->a13 = a13; l->s13 = NaN; + geod_genposition(l, GEOD_ARCMODE, l->a13, 0, 0, 0, &l->s13, 0, 0, 0, 0); +} + +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); +} + void geod_position(const struct geod_geodesicline* l, real s12, real* plat2, real* plon2, real* pazi2) { geod_genposition(l, FALSE, s12, plat2, plon2, pazi2, 0, 0, 0, 0, 0); @@ -630,23 +674,26 @@ void geod_direct(const struct geod_geodesic* g, 0, 0, 0, 0, 0); } -real geod_geninverse(const struct geod_geodesic* g, - real lat1, real lon1, real lat2, real lon2, - real* ps12, real* pazi1, real* pazi2, - real* pm12, real* pM12, real* pM21, real* pS12) { - real s12 = 0, azi1 = 0, azi2 = 0, m12 = 0, M12 = 0, M21 = 0, S12 = 0; - real lon12; +static real geod_geninverse_int(const struct geod_geodesic* g, + real lat1, real lon1, real lat2, real lon2, + real* ps12, + real* psalp1, real* pcalp1, + real* psalp2, real* pcalp2, + real* pm12, real* pM12, real* pM21, + real* pS12) { + real s12 = 0, m12 = 0, M12 = 0, M21 = 0, S12 = 0; + real lon12, lon12s; int latsign, lonsign, swapp; real sbet1, cbet1, sbet2, cbet2, s12x = 0, m12x = 0; real dn1, dn2, lam12, slam12, clam12; real a12 = 0, sig12, calp1 = 0, salp1 = 0, calp2 = 0, salp2 = 0; real Ca[nC]; boolx meridian; - real omg12 = 0; + /* somg12 > 1 marks that it needs to be calculated */ + real omg12 = 0, somg12 = 2, comg12 = 0; unsigned outmask = (ps12 ? GEOD_DISTANCE : 0U) | - (pazi1 || pazi2 ? GEOD_AZIMUTH : 0U) | (pm12 ? GEOD_REDUCEDLENGTH : 0U) | (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) | (pS12 ? GEOD_AREA : 0U); @@ -655,16 +702,25 @@ real geod_geninverse(const struct geod_geodesic* g, /* Compute longitude difference (AngDiff does this carefully). Result is * in [-180, 180] but -180 is only for west-going geodesics. 180 is for * east-going and meridional geodesics. */ - /* If very close to being on the same half-meridian, then make it so. */ - lon12 = AngRound(AngDiff(lon1, lon2)); + lon12 = AngDiff(lon1, lon2, &lon12s); /* Make longitude difference positive. */ lonsign = lon12 >= 0 ? 1 : -1; - lon12 *= lonsign; + /* If very close to being on the same half-meridian, then make it so. */ + lon12 = lonsign * AngRound(lon12); + lon12s = AngRound((180 - lon12) - lonsign * lon12s); + lam12 = lon12 * degree; + if (lon12 > 90) { + sincosdx(lon12s, &slam12, &clam12); + clam12 = -clam12; + } else + sincosdx(lon12, &slam12, &clam12); + /* If really close to the equator, treat as on equator. */ lat1 = AngRound(LatFix(lat1)); lat2 = AngRound(LatFix(lat2)); - /* Swap points so that point with higher (abs) latitude is point 1 */ - swapp = fabs(lat1) >= fabs(lat2) ? 1 : -1; + /* Swap points so that point with higher (abs) latitude is point 1 + * If one latitude is a nan, then it becomes lat1. */ + swapp = fabs(lat1) < fabs(lat2) ? -1 : 1; if (swapp < 0) { lonsign *= -1; swapx(&lat1, &lat2); @@ -712,9 +768,6 @@ real geod_geninverse(const struct geod_geodesic* g, dn1 = sqrt(1 + g->ep2 * sq(sbet1)); dn2 = sqrt(1 + g->ep2 * sq(sbet2)); - lam12 = lon12 * degree; - sincosdx(lon12, &slam12, &clam12); - meridian = lat1 == -90 || slam12 == 0; if (meridian) { @@ -731,7 +784,7 @@ real geod_geninverse(const struct geod_geodesic* g, ssig2 = sbet2; csig2 = calp2 * cbet2; /* sig12 = sig2 - sig1 */ - sig12 = atan2(maxx(csig1 * ssig2 - ssig1 * csig2, (real)(0)), + sig12 = atan2(maxx((real)(0), csig1 * ssig2 - ssig1 * csig2), csig1 * csig2 + ssig1 * ssig2); Lengths(g, g->n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2, &s12x, &m12x, 0, @@ -760,7 +813,7 @@ real geod_geninverse(const struct geod_geodesic* g, if (!meridian && sbet1 == 0 && /* and sbet2 == 0 */ /* Mimic the way Lambda12 works with calp1 = 0 */ - (g->f <= 0 || lam12 <= pi - g->f * pi)) { + (g->f <= 0 || lon12s >= g->f * 180)) { /* Geodesic runs along equator */ calp1 = calp2 = 0; salp1 = salp2 = 1; @@ -779,7 +832,7 @@ real geod_geninverse(const struct geod_geodesic* g, /* Figure a starting point for Newton's method */ real dnm = 0; sig12 = InverseStart(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, - lam12, + lam12, slam12, clam12, &salp1, &calp1, &salp2, &calp2, &dnm, Ca); @@ -813,13 +866,13 @@ real geod_geninverse(const struct geod_geodesic* g, /* the WGS84 test set: mean = 1.47, sd = 1.25, max = 16 * WGS84 and random input: mean = 2.85, sd = 0.60 */ real dv = 0, - v = (Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1, + v = Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1, + slam12, clam12, &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2, - &eps, &omg12, numit < maxit1, &dv, Ca) - - lam12); + &eps, &somg12, &comg12, 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 : 2) * tol0)) break; + if (tripb || !(fabs(v) >= (tripn ? 8 : 1) * tol0)) break; /* Update bracketing values */ if (v > 0 && (numit > maxit1 || calp1/salp1 > calp1b/salp1b)) { salp1b = salp1; calp1b = calp1; } @@ -864,7 +917,6 @@ real geod_geninverse(const struct geod_geodesic* g, m12x *= g->b; s12x *= g->b; a12 = sig12 / degree; - omg12 = lam12 - omg12; } } @@ -900,15 +952,22 @@ real geod_geninverse(const struct geod_geodesic* g, /* Avoid problems with indeterminate sig1, sig2 on equator */ S12 = 0; + if (!meridian) { + if (somg12 > 1) { + somg12 = sin(omg12); comg12 = cos(omg12); + } else + norm2(&somg12, &comg12); + } + if (!meridian && - omg12 < (real)(0.75) * pi && /* Long difference too big */ - sbet2 - sbet1 < (real)(1.75)) { /* Lat difference too big */ + /* omg12 < 3/4 * pi */ + comg12 > -(real)(0.7071) && /* Long difference not too big */ + sbet2 - sbet1 < (real)(1.75)) { /* Lat difference not too big */ /* Use tan(Gamma/2) = tan(omg12/2) * * (tan(bet1/2)+tan(bet2/2))/(1+tan(bet1/2)*tan(bet2/2)) * with tan(x/2) = sin(x)/(1+cos(x)) */ real - somg12 = sin(omg12), domg12 = 1 + cos(omg12), - dbet1 = 1 + cbet1, dbet2 = 1 + cbet2; + domg12 = 1 + comg12, dbet1 = 1 + cbet1, dbet2 = 1 + cbet2; alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ), domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) ); } else { @@ -943,18 +1002,13 @@ real geod_geninverse(const struct geod_geodesic* g, salp1 *= swapp * lonsign; calp1 *= swapp * latsign; salp2 *= swapp * lonsign; calp2 *= swapp * latsign; - if (outmask & GEOD_AZIMUTH) { - /* minus signs give range [-180, 180). 0- converts -0 to +0. */ - azi1 = atan2dx(salp1, calp1); - azi2 = atan2dx(salp2, calp2); - } + if (psalp1) *psalp1 = salp1; + if (pcalp1) *pcalp1 = calp1; + if (psalp2) *psalp2 = salp2; + if (pcalp2) *pcalp2 = calp2; if (outmask & GEOD_DISTANCE) *ps12 = s12; - if (outmask & GEOD_AZIMUTH) { - if (pazi1) *pazi1 = azi1; - if (pazi2) *pazi2 = azi2; - } if (outmask & GEOD_REDUCEDLENGTH) *pm12 = m12; if (outmask & GEOD_GEODESICSCALE) { @@ -968,6 +1022,35 @@ real geod_geninverse(const struct geod_geodesic* g, return a12; } +real geod_geninverse(const struct geod_geodesic* g, + real lat1, real lon1, real lat2, real lon2, + real* ps12, real* pazi1, real* pazi2, + real* pm12, real* pM12, real* pM21, real* pS12) { + real salp1, calp1, salp2, calp2, + a12 = geod_geninverse_int(g, lat1, lon1, lat2, lon2, ps12, + &salp1, &calp1, &salp2, &calp2, + pm12, pM12, pM21, pS12); + if (pazi1) *pazi1 = atan2dx(salp1, calp1); + if (pazi2) *pazi2 = atan2dx(salp2, calp2); + return a12; +} + +void geod_inverseline(struct geod_geodesicline* l, + const struct geod_geodesic* g, + real lat1, real lon1, real lat2, real lon2, + unsigned caps) { + real salp1, calp1, + a12 = geod_geninverse_int(g, lat1, lon1, lat2, lon2, 0, + &salp1, &calp1, 0, 0, + 0, 0, 0, 0), + azi1 = atan2dx(salp1, calp1); + caps = caps ? caps : GEOD_DISTANCE_IN | GEOD_LONGITUDE; + /* Ensure that a12 can be converted to a distance */ + if (caps & (OUT_ALL & GEOD_DISTANCE_IN)) caps |= GEOD_DISTANCE; + geod_lineinit_int(l, g, lat1, lon1, azi1, salp1, calp1, caps); + geod_setarc(l, a12); +} + void geod_inverse(const struct geod_geodesic* g, real lat1, real lon1, real lat2, real lon2, real* ps12, real* pazi1, real* pazi2) { @@ -1112,7 +1195,7 @@ real Astroid(real x, real y) { real InverseStart(const struct geod_geodesic* g, real sbet1, real cbet1, real dn1, real sbet2, real cbet2, real dn2, - real lam12, + real lam12, real slam12, real clam12, real* psalp1, real* pcalp1, /* Only updated if return val >= 0 */ real* psalp2, real* pcalp2, @@ -1133,7 +1216,7 @@ real InverseStart(const struct geod_geodesic* g, real sbet12a; boolx shortline = cbet12 >= 0 && sbet12 < (real)(0.5) && cbet2 * lam12 < (real)(0.5); - real omg12 = lam12, somg12, comg12, ssig12, csig12; + real somg12, comg12, ssig12, csig12; #if defined(__GNUC__) && __GNUC__ == 4 && \ (__GNUC_MINOR__ < 6 || defined(__MINGW32__)) /* Volatile declaration needed to fix inverse cases @@ -1151,14 +1234,16 @@ real InverseStart(const struct geod_geodesic* g, sbet12a = sbet2 * cbet1 + cbet2 * sbet1; #endif if (shortline) { - real sbetm2 = sq(sbet1 + sbet2); + real sbetm2 = sq(sbet1 + sbet2), omg12; /* sin((bet1+bet2)/2)^2 * = (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2) */ sbetm2 /= sbetm2 + sq(cbet1 + cbet2); dnm = sqrt(1 + g->ep2 * sbetm2); - omg12 /= g->f1 * dnm; + omg12 = lam12 / (g->f1 * dnm); + somg12 = sin(omg12); comg12 = cos(omg12); + } else { + somg12 = slam12; comg12 = clam12; } - somg12 = sin(omg12); comg12 = cos(omg12); salp1 = cbet2 * somg12; calp1 = comg12 >= 0 ? @@ -1188,6 +1273,7 @@ real InverseStart(const struct geod_geodesic* g, * 56.320923501171 0 -56.320923501171 179.664747671772880215 * which otherwise fails with g++ 4.4.4 x86 -O3 */ volatile real x; + real lam12x = atan2(-slam12, -clam12); /* lam12 - pi */ if (g->f >= 0) { /* In fact f == 0 does not get here */ /* x = dlong, y = dlat */ { @@ -1198,7 +1284,7 @@ real InverseStart(const struct geod_geodesic* g, } betscale = lamscale * cbet1; - x = (lam12 - pi) / lamscale; + x = lam12x / lamscale; y = sbet12a / betscale; } else { /* f < 0 */ /* x = dlat, y = dlong */ @@ -1215,7 +1301,7 @@ real InverseStart(const struct geod_geodesic* g, betscale = x < -(real)(0.01) ? sbet12a / x : -g->f * sq(cbet1) * pi; lamscale = betscale / cbet1; - y = (lam12 - pi) / lamscale; + y = lam12x / lamscale; } if (y > -tol1 && x > -1 - xthresh) { @@ -1292,19 +1378,22 @@ real Lambda12(const struct geod_geodesic* g, real sbet1, real cbet1, real dn1, real sbet2, real cbet2, real dn2, real salp1, real calp1, + real slam120, real clam120, real* psalp2, real* pcalp2, real* psig12, real* pssig1, real* pcsig1, real* pssig2, real* pcsig2, - real* peps, real* pdomg12, + real* peps, + real* psomg12, real* pcomg12, boolx diffp, real* pdlam12, /* Scratch area of the right size */ real Ca[]) { real salp2 = 0, calp2 = 0, sig12 = 0, - ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0, domg12 = 0, dlam12 = 0; + ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0, + somg12 = 0, comg12 = 0, dlam12 = 0; real salp0, calp0; - real somg1, comg1, somg2, comg2, omg12, lam12; - real B312, h0, k2; + real somg1, comg1, somg2, comg2, lam12; + real B312, eta, k2; if (sbet1 == 0 && calp1 == 0) /* Break degeneracy of equatorial line. This case has already been @@ -1345,20 +1434,21 @@ real Lambda12(const struct geod_geodesic* g, /* norm2(&somg2, &comg2); -- don't need to normalize! */ /* sig12 = sig2 - sig1, limit to [0, pi] */ - sig12 = atan2(maxx(csig1 * ssig2 - ssig1 * csig2, (real)(0)), + sig12 = atan2(maxx((real)(0), csig1 * ssig2 - ssig1 * csig2), csig1 * csig2 + ssig1 * ssig2); /* omg12 = omg2 - omg1, limit to [0, pi] */ - omg12 = atan2(maxx(comg1 * somg2 - somg1 * comg2, (real)(0)), - comg1 * comg2 + somg1 * somg2); + somg12 = maxx((real)(0), comg1 * somg2 - somg1 * comg2); + comg12 = comg1 * comg2 + somg1 * somg2; + /* eta = omg12 - lam120 */ + eta = atan2(somg12 * clam120 - comg12 * slam120, + comg12 * clam120 + somg12 * slam120); k2 = sq(calp0) * g->ep2; eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2); C3f(g, eps, Ca); B312 = (SinCosSeries(TRUE, ssig2, csig2, Ca, nC3-1) - SinCosSeries(TRUE, ssig1, csig1, Ca, nC3-1)); - h0 = -g->f * A3f(g, eps); - domg12 = salp0 * h0 * (sig12 + B312); - lam12 = omg12 + domg12; + lam12 = eta - g->f * A3f(g, eps) * salp0 * (sig12 + B312); if (diffp) { if (calp2 == 0) @@ -1378,7 +1468,8 @@ real Lambda12(const struct geod_geodesic* g, *pssig2 = ssig2; *pcsig2 = csig2; *peps = eps; - *pdomg12 = domg12; + *psomg12 = somg12; + *pcomg12 = comg12; if (diffp) *pdlam12 = dlam12; @@ -1653,7 +1744,7 @@ int transit(real lon1, real lon2) { /* Compute lon12 the same way as Geodesic::Inverse. */ lon1 = AngNormalize(lon1); lon2 = AngNormalize(lon2); - lon12 = AngDiff(lon1, lon2); + lon12 = AngDiff(lon1, lon2, 0); return lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 : (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0); } @@ -1699,8 +1790,12 @@ void accneg(real s[]) { } void geod_polygon_init(struct geod_polygon* p, boolx polylinep) { - p->lat0 = p->lon0 = p->lat = p->lon = NaN; p->polyline = (polylinep != 0); + geod_polygon_clear(p); +} + +void geod_polygon_clear(struct geod_polygon* p) { + p->lat0 = p->lon0 = p->lat = p->lon = NaN; accini(p->P); accini(p->A); p->num = p->crossings = 0; diff --git a/src/geodesic.h b/src/geodesic.h index 7bd8270f..80ffac9b 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -1,6 +1,6 @@ /** * \file geodesic.h - * \brief Header for the geodesic routines in C + * \brief API for the geodesic routines in C * * This an implementation in C of the geodesic algorithms described in * - C. F. F. Karney, @@ -9,7 +9,7 @@ * J. Geodesy <b>87</b>, 43--55 (2013); * DOI: <a href="https://dx.doi.org/10.1007/s00190-012-0578-z"> * 10.1007/s00190-012-0578-z</a>; - * addenda: <a href="http://geographiclib.sf.net/geod-addenda.html"> + * addenda: <a href="http://geographiclib.sourceforge.net/geod-addenda.html"> * geod-addenda.html</a>. * . * The principal advantages of these algorithms over previous ones (e.g., @@ -75,30 +75,29 @@ * (obviously) uniquely defined. However, in a few special cases there are * multiple azimuths which yield the same shortest distance. Here is a * catalog of those cases: - * - \e lat1 = −\e lat2 (with neither point at a pole). If \e azi1 = - * \e azi2, the geodesic is unique. Otherwise there are two geodesics - * and the second one is obtained by setting [\e azi1, \e azi2] = [\e - * azi2, \e azi1], [\e M12, \e M21] = [\e M21, \e M12], \e S12 = - * −\e S12. (This occurs when the longitude difference is near - * ±180° for oblate ellipsoids.) - * - \e lon2 = \e lon1 ± 180° (with neither point at a pole). - * If \e azi1 = 0° or ±180°, the geodesic is unique. - * Otherwise there are two geodesics and the second one is obtained by - * setting [\e azi1, \e azi2] = [−\e azi1, −\e azi2], \e S12 - * = −\e S12. (This occurs when \e lat2 is near −\e lat1 for + * - \e lat1 = −\e lat2 (with neither point at a pole). If \e azi1 = \e + * azi2, the geodesic is unique. Otherwise there are two geodesics and the + * second one is obtained by setting [\e azi1, \e azi2] → [\e azi2, \e + * azi1], [\e M12, \e M21] → [\e M21, \e M12], \e S12 → −\e + * S12. (This occurs when the longitude difference is near ±180° + * for oblate ellipsoids.) + * - \e lon2 = \e lon1 ± 180° (with neither point at a pole). If \e + * azi1 = 0° or ±180°, the geodesic is unique. Otherwise + * there are two geodesics and the second one is obtained by setting [\e + * azi1, \e azi2] → [−\e azi1, −\e azi2], \e S12 → + * −\e S12. (This occurs when \e lat2 is near −\e lat1 for * prolate ellipsoids.) - * - Points 1 and 2 at opposite poles. There are infinitely many - * geodesics which can be generated by setting [\e azi1, \e azi2] = - * [\e azi1, \e azi2] + [\e d, −\e d], for arbitrary \e d. (For - * spheres, this prescription applies when points 1 and 2 are - * antipodal.) - * - \e s12 = 0 (coincident points). There are infinitely many geodesics - * which can be generated by setting [\e azi1, \e azi2] = [\e azi1, \e - * azi2] + [\e d, \e d], for arbitrary \e d. + * - Points 1 and 2 at opposite poles. There are infinitely many geodesics + * which can be generated by setting [\e azi1, \e azi2] → [\e azi1, \e + * azi2] + [\e d, −\e d], for arbitrary \e d. (For spheres, this + * prescription applies when points 1 and 2 are antipodal.) + * - \e s12 = 0 (coincident points). There are infinitely many geodesics which + * can be generated by setting [\e azi1, \e azi2] → [\e azi1, \e azi2] + + * [\e d, \e d], for arbitrary \e d. * * These routines are a simple transcription of the corresponding C++ classes - * in <a href="http://geographiclib.sf.net"> GeographicLib</a>. The "class - * data" is represented by the structs geod_geodesic, geod_geodesicline, + * in <a href="http://geographiclib.sourceforge.net"> GeographicLib</a>. The + * "class data" is represented by the structs geod_geodesic, geod_geodesicline, * geod_polygon and pointers to these objects are passed as initial arguments * to the member functions. Most of the internal comments have been retained. * However, in the process of transcription some documentation has been lost @@ -108,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-2015) <charles@karney.com> and licensed + * Copyright (c) Charles Karney (2012-2016) <charles@karney.com> and licensed * under the MIT/X11 License. For more information, see * http://geographiclib.sourceforge.net/ * * This library was distributed with - * <a href="../index.html">GeographicLib</a> 1.44. + * <a href="../index.html">GeographicLib</a> 1.46. **********************************************************************/ #if !defined(GEODESIC_H) @@ -128,7 +127,7 @@ * The minor version of the geodesic library. (This tracks the version of * GeographicLib.) **********************************************************************/ -#define GEODESIC_VERSION_MINOR 44 +#define GEODESIC_VERSION_MINOR 46 /** * The patch level of the geodesic library. (This tracks the version of * GeographicLib.) @@ -177,7 +176,8 @@ extern "C" { /** * The struct containing information about a single geodesic. This must be - * initialized by geod_lineinit() before use. + * initialized by geod_lineinit(), geod_directline(), geod_gendirectline(), + * or geod_inverseline() before use. **********************************************************************/ struct geod_geodesicline { double lat1; /**< the starting latitude */ @@ -185,9 +185,13 @@ extern "C" { double azi1; /**< the starting azimuth */ double a; /**< the equatorial radius */ double f; /**< the flattening */ + double salp1; /**< sine of \e azi1 */ + double calp1; /**< cosine of \e azi1 */ + double a13; /**< arc length to reference point */ + double s13; /**< distance to reference point */ /**< @cond SKIP */ double b, c2, f1, salp0, calp0, k2, - salp1, calp1, ssig1, csig1, dn1, stau1, ctau1, somg1, comg1, + ssig1, csig1, dn1, stau1, ctau1, somg1, comg1, A1m1, A2m1, A3c, B11, B21, B31, A4, B41; double C1a[6+1], C1pa[6+1], C2a[6+1], C3a[6], C4a[6]; /**< @endcond */ @@ -223,46 +227,6 @@ extern "C" { void geod_init(struct geod_geodesic* g, double a, double f); /** - * Initialize a geod_geodesicline object. - * - * @param[out] l a pointer to the object to be initialized. - * @param[in] g a pointer to the geod_geodesic object specifying the - * ellipsoid. - * @param[in] lat1 latitude of point 1 (degrees). - * @param[in] lon1 longitude of point 1 (degrees). - * @param[in] azi1 azimuth at point 1 (degrees). - * @param[in] caps bitor'ed combination of geod_mask() values specifying the - * capabilities the geod_geodesicline object should possess, i.e., which - * quantities can be returned in calls to geod_position() and - * geod_genposition(). - * - * \e g must have been initialized with a call to geod_init(). \e lat1 - * should be in the range [−90°, 90°]. - * - * The geod_mask values are [see geod_mask()]: - * - \e caps |= GEOD_LATITUDE for the latitude \e lat2; this is - * added automatically, - * - \e caps |= GEOD_LONGITUDE for the latitude \e lon2, - * - \e caps |= GEOD_AZIMUTH for the latitude \e azi2; this is - * added automatically, - * - \e caps |= GEOD_DISTANCE for the distance \e s12, - * - \e caps |= GEOD_REDUCEDLENGTH for the reduced length \e m12, - * - \e caps |= GEOD_GEODESICSCALE for the geodesic scales \e M12 - * and \e M21, - * - \e caps |= GEOD_AREA for the area \e S12, - * - \e caps |= GEOD_DISTANCE_IN permits the length of the - * geodesic to be given in terms of \e s12; without this capability the - * length can only be specified in terms of arc length. - * . - * A value of \e caps = 0 is treated as GEOD_LATITUDE | GEOD_LONGITUDE | - * GEOD_AZIMUTH | GEOD_DISTANCE_IN (to support the solution of the "standard" - * direct problem). - **********************************************************************/ - void geod_lineinit(struct geod_geodesicline* l, - const struct geod_geodesic* g, - double lat1, double lon1, double azi1, unsigned caps); - - /** * Solve the direct geodesic problem. * * @param[in] g a pointer to the geod_geodesic object specifying the @@ -270,7 +234,7 @@ extern "C" { * @param[in] lat1 latitude of point 1 (degrees). * @param[in] lon1 longitude of point 1 (degrees). * @param[in] azi1 azimuth at point 1 (degrees). - * @param[in] s12 distance between point 1 and point 2 (meters); it can be + * @param[in] s12 distance from point 1 to point 2 (meters); it can be * negative. * @param[out] plat2 pointer to the latitude of point 2 (degrees). * @param[out] plon2 pointer to the longitude of point 2 (degrees). @@ -303,6 +267,51 @@ extern "C" { double* plat2, double* plon2, double* pazi2); /** + * The general direct geodesic problem. + * + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] lat1 latitude of point 1 (degrees). + * @param[in] lon1 longitude of point 1 (degrees). + * @param[in] azi1 azimuth at point 1 (degrees). + * @param[in] flags bitor'ed combination of geod_flags(); \e flags & + * GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags & + * GEOD_LONG_UNROLL "unrolls" \e lon2. + * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the distance + * from point 1 to point 2 (meters); otherwise it is the arc length + * from point 1 to point 2 (degrees); it can be negative. + * @param[out] plat2 pointer to the latitude of point 2 (degrees). + * @param[out] plon2 pointer to the longitude of point 2 (degrees). + * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). + * @param[out] ps12 pointer to the distance from point 1 to point 2 + * (meters). + * @param[out] pm12 pointer to the reduced length of geodesic (meters). + * @param[out] pM12 pointer to the geodesic scale of point 2 relative to + * point 1 (dimensionless). + * @param[out] pM21 pointer to the geodesic scale of point 1 relative to + * point 2 (dimensionless). + * @param[out] pS12 pointer to the area under the geodesic + * (meters<sup>2</sup>). + * @return \e a12 arc length from point 1 to point 2 (degrees). + * + * \e g must have been initialized with a call to geod_init(). \e lat1 + * should be in the range [−90°, 90°]. The function value \e + * a12 equals \e s12_a12 if \e flags & GEOD_ARCMODE. Any of the "return" + * arguments, \e plat2, etc., may be replaced by 0, if you do not need some + * quantities computed. + * + * With \e flags & GEOD_LONG_UNROLL bit set, the longitude is "unrolled" so + * that the quantity \e lon2 − \e lon1 indicates how many times and in + * what sense the geodesic encircles the ellipsoid. + **********************************************************************/ + double 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); + + /** * Solve the inverse geodesic problem. * * @param[in] g a pointer to the geod_geodesic object specifying the @@ -311,7 +320,7 @@ extern "C" { * @param[in] lon1 longitude of point 1 (degrees). * @param[in] lat2 latitude of point 2 (degrees). * @param[in] lon2 longitude of point 2 (degrees). - * @param[out] ps12 pointer to the distance between point 1 and point 2 + * @param[out] ps12 pointer to the distance from point 1 to point 2 * (meters). * @param[out] pazi1 pointer to the azimuth at point 1 (degrees). * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). @@ -345,22 +354,180 @@ extern "C" { double* ps12, double* pazi1, double* pazi2); /** + * The general inverse geodesic calculation. + * + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] lat1 latitude of point 1 (degrees). + * @param[in] lon1 longitude of point 1 (degrees). + * @param[in] lat2 latitude of point 2 (degrees). + * @param[in] lon2 longitude of point 2 (degrees). + * @param[out] ps12 pointer to the distance from point 1 to point 2 + * (meters). + * @param[out] pazi1 pointer to the azimuth at point 1 (degrees). + * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). + * @param[out] pm12 pointer to the reduced length of geodesic (meters). + * @param[out] pM12 pointer to the geodesic scale of point 2 relative to + * point 1 (dimensionless). + * @param[out] pM21 pointer to the geodesic scale of point 1 relative to + * point 2 (dimensionless). + * @param[out] pS12 pointer to the area under the geodesic + * (meters<sup>2</sup>). + * @return \e a12 arc length from point 1 to point 2 (degrees). + * + * \e g must have been initialized with a call to geod_init(). \e lat1 and + * \e lat2 should be in the range [−90°, 90°]. Any of the + * "return" arguments \e ps12, etc., may be replaced by 0, if you do not need + * some quantities computed. + **********************************************************************/ + double 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); + + /** + * Initialize a geod_geodesicline object. + * + * @param[out] l a pointer to the object to be initialized. + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] lat1 latitude of point 1 (degrees). + * @param[in] lon1 longitude of point 1 (degrees). + * @param[in] azi1 azimuth at point 1 (degrees). + * @param[in] caps bitor'ed combination of geod_mask() values specifying the + * capabilities the geod_geodesicline object should possess, i.e., which + * quantities can be returned in calls to geod_position() and + * geod_genposition(). + * + * \e g must have been initialized with a call to geod_init(). \e lat1 + * should be in the range [−90°, 90°]. + * + * The geod_mask values are [see geod_mask()]: + * - \e caps |= GEOD_LATITUDE for the latitude \e lat2; this is + * added automatically, + * - \e caps |= GEOD_LONGITUDE for the latitude \e lon2, + * - \e caps |= GEOD_AZIMUTH for the latitude \e azi2; this is + * added automatically, + * - \e caps |= GEOD_DISTANCE for the distance \e s12, + * - \e caps |= GEOD_REDUCEDLENGTH for the reduced length \e m12, + * - \e caps |= GEOD_GEODESICSCALE for the geodesic scales \e M12 + * and \e M21, + * - \e caps |= GEOD_AREA for the area \e S12, + * - \e caps |= GEOD_DISTANCE_IN permits the length of the + * geodesic to be given in terms of \e s12; without this capability the + * length can only be specified in terms of arc length. + * . + * A value of \e caps = 0 is treated as GEOD_LATITUDE | GEOD_LONGITUDE | + * GEOD_AZIMUTH | GEOD_DISTANCE_IN (to support the solution of the "standard" + * direct problem). + * + * When initialized by this function, point 3 is undefined (l->s13 = l->a13 = + * NaN). + **********************************************************************/ + void geod_lineinit(struct geod_geodesicline* l, + const struct geod_geodesic* g, + double lat1, double lon1, double azi1, unsigned caps); + + /** + * Initialize a geod_geodesicline object in terms of the direct geodesic + * problem. + * + * @param[out] l a pointer to the object to be initialized. + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] lat1 latitude of point 1 (degrees). + * @param[in] lon1 longitude of point 1 (degrees). + * @param[in] azi1 azimuth at point 1 (degrees). + * @param[in] s12 distance from point 1 to point 2 (meters); it can be + * negative. + * @param[in] caps bitor'ed combination of geod_mask() values specifying the + * capabilities the geod_geodesicline object should possess, i.e., which + * quantities can be returned in calls to geod_position() and + * geod_genposition(). + * + * This function sets point 3 of the geod_geodesicline to correspond to point + * 2 of the direct geodesic problem. See geod_lineinit() for more + * informaion. + **********************************************************************/ + void geod_directline(struct geod_geodesicline* l, + 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 + * problem spacified in terms of either distance or arc length. + * + * @param[out] l a pointer to the object to be initialized. + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] lat1 latitude of point 1 (degrees). + * @param[in] lon1 longitude of point 1 (degrees). + * @param[in] azi1 azimuth at point 1 (degrees). + * @param[in] flags either GEOD_NOFLAGS or GEOD_ARCMODE to determining the + * meaning of the \e s12_a12. + * @param[in] s12_a12 if \e flags = GEOD_NOFLAGS, this is the distance + * from point 1 to point 2 (meters); if \e flags = GEOD_ARCMODE, it is + * the arc length from point 1 to point 2 (degrees); it can be + * negative. + * @param[in] caps bitor'ed combination of geod_mask() values specifying the + * capabilities the geod_geodesicline object should possess, i.e., which + * quantities can be returned in calls to geod_position() and + * geod_genposition(). + * + * This function sets point 3 of the geod_geodesicline to correspond to point + * 2 of the direct geodesic problem. See geod_lineinit() for more + * informaion. + **********************************************************************/ + void geod_gendirectline(struct geod_geodesicline* l, + 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 + * problem. + * + * @param[out] l a pointer to the object to be initialized. + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] lat1 latitude of point 1 (degrees). + * @param[in] lon1 longitude of point 1 (degrees). + * @param[in] lat2 latitude of point 2 (degrees). + * @param[in] lon2 longitude of point 2 (degrees). + * @param[in] caps bitor'ed combination of geod_mask() values specifying the + * capabilities the geod_geodesicline object should possess, i.e., which + * quantities can be returned in calls to geod_position() and + * geod_genposition(). + * + * This function sets point 3 of the geod_geodesicline to correspond to point + * 2 of the inverse geodesic problem. See geod_lineinit() for more + * informaion. + **********************************************************************/ + void geod_inverseline(struct geod_geodesicline* l, + const struct geod_geodesic* g, + double lat1, double lon1, double lat2, double lon2, + unsigned caps); + + /** * Compute the position along a geod_geodesicline. * * @param[in] l a pointer to the geod_geodesicline object specifying the * geodesic line. - * @param[in] s12 distance between point 1 and point 2 (meters); it can be + * @param[in] s12 distance from point 1 to point 2 (meters); it can be * negative. * @param[out] plat2 pointer to the latitude of point 2 (degrees). * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires * that \e l was initialized with \e caps |= GEOD_LONGITUDE. * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). * - * \e l must have been initialized with a call 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 need some quantities - * computed. + * \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 + * need some quantities computed. * * Example, compute way points between JFK and Singapore Changi Airport * the "obvious" way using geod_direct(): @@ -379,13 +546,12 @@ extern "C" { @code{.c} struct geod_geodesic g; struct geod_geodesicline l; - double s12, azi1, lat[101],lon[101]; + double lat[101],lon[101]; int i; geod_init(&g, 6378137, 1/298.257223563); - geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0); - geod_lineinit(&l, &g, 40.64, -73.78, azi1, 0); - for (i = 0; i < 101; ++i) { - geod_position(&l, i * s12 * 0.01, lat + i, lon + i, 0); + geod_inverseline(&l, &g, 40.64, -73.78, 1.36, 103.99, 0); + for (i = 0; i <= 100; ++i) { + geod_position(&l, i * l.s13 * 0.01, lat + i, lon + i, 0); printf("%.5f %.5f\n", lat[i], lon[i]); } @endcode @@ -394,84 +560,6 @@ extern "C" { double* plat2, double* plon2, double* pazi2); /** - * The general direct geodesic problem. - * - * @param[in] g a pointer to the geod_geodesic object specifying the - * ellipsoid. - * @param[in] lat1 latitude of point 1 (degrees). - * @param[in] lon1 longitude of point 1 (degrees). - * @param[in] azi1 azimuth at point 1 (degrees). - * @param[in] flags bitor'ed combination of geod_flags(); \e flags & - * GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags & - * GEOD_LONG_UNROLL "unrolls" \e lon2. - * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the distance - * between point 1 and point 2 (meters); otherwise it is the arc length - * between point 1 and point 2 (degrees); it can be negative. - * @param[out] plat2 pointer to the latitude of point 2 (degrees). - * @param[out] plon2 pointer to the longitude of point 2 (degrees). - * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). - * @param[out] ps12 pointer to the distance between point 1 and point 2 - * (meters). - * @param[out] pm12 pointer to the reduced length of geodesic (meters). - * @param[out] pM12 pointer to the geodesic scale of point 2 relative to - * point 1 (dimensionless). - * @param[out] pM21 pointer to the geodesic scale of point 1 relative to - * point 2 (dimensionless). - * @param[out] pS12 pointer to the area under the geodesic - * (meters<sup>2</sup>). - * @return \e a12 arc length of between point 1 and point 2 (degrees). - * - * \e g must have been initialized with a call to geod_init(). \e lat1 - * should be in the range [−90°, 90°]. The function value \e - * a12 equals \e s12_a12 if \e flags & GEOD_ARCMODE. Any of the "return" - * arguments, \e plat2, etc., may be replaced by 0, if you do not need some - * quantities computed. - * - * With \e flags & GEOD_LONG_UNROLL bit set, the longitude is "unrolled" so - * that the quantity \e lon2 − \e lon1 indicates how many times and in - * what sense the geodesic encircles the ellipsoid. - **********************************************************************/ - double 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); - - /** - * The general inverse geodesic calculation. - * - * @param[in] g a pointer to the geod_geodesic object specifying the - * ellipsoid. - * @param[in] lat1 latitude of point 1 (degrees). - * @param[in] lon1 longitude of point 1 (degrees). - * @param[in] lat2 latitude of point 2 (degrees). - * @param[in] lon2 longitude of point 2 (degrees). - * @param[out] ps12 pointer to the distance between point 1 and point 2 - * (meters). - * @param[out] pazi1 pointer to the azimuth at point 1 (degrees). - * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). - * @param[out] pm12 pointer to the reduced length of geodesic (meters). - * @param[out] pM12 pointer to the geodesic scale of point 2 relative to - * point 1 (dimensionless). - * @param[out] pM21 pointer to the geodesic scale of point 1 relative to - * point 2 (dimensionless). - * @param[out] pS12 pointer to the area under the geodesic - * (meters<sup>2</sup>). - * @return \e a12 arc length of between point 1 and point 2 (degrees). - * - * \e g must have been initialized with a call to geod_init(). \e lat1 and - * \e lat2 should be in the range [−90°, 90°]. Any of the - * "return" arguments \e ps12, etc., may be replaced by 0, if you do not need - * some quantities computed. - **********************************************************************/ - double 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); - - /** * The general position function. * * @param[in] l a pointer to the geod_geodesicline object specifying the @@ -481,14 +569,14 @@ extern "C" { * GEOD_LONG_UNROLL "unrolls" \e lon2; if \e flags & GEOD_ARCMODE is 0, * then \e l must have been initialized with \e caps |= GEOD_DISTANCE_IN. * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the - * distance between point 1 and point 2 (meters); otherwise it is the - * arc length between point 1 and point 2 (degrees); it can be + * distance from point 1 to point 2 (meters); otherwise it is the + * arc length from point 1 to point 2 (degrees); it can be * negative. * @param[out] plat2 pointer to the latitude of point 2 (degrees). * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires * that \e l was initialized with \e caps |= GEOD_LONGITUDE. * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). - * @param[out] ps12 pointer to the distance between point 1 and point 2 + * @param[out] ps12 pointer to the distance from point 1 to point 2 * (meters); requires that \e l was initialized with \e caps |= * GEOD_DISTANCE. * @param[out] pm12 pointer to the reduced length of geodesic (meters); @@ -502,7 +590,7 @@ extern "C" { * @param[out] pS12 pointer to the area under the geodesic * (meters<sup>2</sup>); requires that \e l was initialized with \e caps |= * GEOD_AREA. - * @return \e a12 arc length of between point 1 and point 2 (degrees). + * @return \e a12 arc length from point 1 to point 2 (degrees). * * \e l must have been initialized with a call to geod_lineinit() with \e * caps |= GEOD_DISTANCE_IN. The value \e azi2 returned is in the range @@ -515,22 +603,21 @@ extern "C" { * that the quantity \e lon2 − \e lon1 indicates how many times and in * what sense the geodesic encircles the ellipsoid. * - * Example, compute way points between JFK and Singapore Changi Airport - * using geod_genposition(). In this example, the points are evenly space in - * arc length (and so only approximately equally space in distance). This is - * faster than using geod_position() would be appropriate if drawing the path - * on a map. + * Example, compute way points between JFK and Singapore Changi Airport using + * geod_genposition(). In this example, the points are evenly space in arc + * length (and so only approximately equally spaced in distance). This is + * faster than using geod_position() and would be appropriate if drawing the + * path on a map. @code{.c} struct geod_geodesic g; struct geod_geodesicline l; - double a12, azi1, lat[101], lon[101]; + double lat[101], lon[101]; int i; geod_init(&g, 6378137, 1/298.257223563); - a12 = geod_geninverse(&g, 40.64, -73.78, 1.36, 103.99, - 0, &azi1, 0, 0, 0, 0, 0); - geod_lineinit(&l, &g, 40.64, -73.78, azi1, GEOD_LATITUDE | GEOD_LONGITUDE); - for (i = 0; i < 101; ++i) { - geod_genposition(&l, 1, i * a12 * 0.01, + geod_inverseline(&l, &g, 40.64, -73.78, 1.36, 103.99, + GEOD_LATITUDE | GEOD_LONGITUDE); + for (i = 0; i <= 100; ++i) { + geod_genposition(&l, GEOD_ARCMODE, i * l.a13 * 0.01, lat + i, lon + i, 0, 0, 0, 0, 0, 0); printf("%.5f %.5f\n", lat[i], lon[i]); } @@ -544,6 +631,36 @@ extern "C" { double* pS12); /** + * Specify position of point 3 in terms of distance. + * + * @param[inout] l a pointer to the geod_geodesicline object. + * @param[in] s13 the distance from point 1 to point 3 (meters); it + * can be negative. + * + * This is only useful if the geod_geodesicline object has been constructed + * with \e caps |= GEOD_DISTANCE_IN. + **********************************************************************/ + void geod_setdistance(struct geod_geodesicline* l, double s13); + + /** + * Specify position of point 3 in terms of either distance or arc length. + * + * @param[inout] l a pointer to the geod_geodesicline object. + * @param[in] flags either GEOD_NOFLAGS or GEOD_ARCMODE to determining the + * meaning of the \e s13_a13. + * @param[in] s13_a13 if \e flags = GEOD_NOFLAGS, this is the distance + * from point 1 to point 3 (meters); if \e flags = GEOD_ARCMODE, it is + * the arc length from point 1 to point 3 (degrees); it can be + * negative. + * + * If flags = GEOD_NOFLAGS, this calls geod_setdistance(). If flags = + * GEOD_ARCMODE, the \e s13 is only set if the geod_geodesicline object has + * been constructed with \e caps |= GEOD_DISTANCE. + **********************************************************************/ + void geod_gensetdistance(struct geod_geodesicline* l, + unsigned flags, double s13_a13); + + /** * Initialize a geod_polygon object. * * @param[out] p a pointer to the object to be initialized. @@ -565,6 +682,13 @@ extern "C" { void geod_polygon_init(struct geod_polygon* p, int polylinep); /** + * Clear the polygon, allowing a new polygon to be started. + * + * @param[in,out] p a pointer to the object to be cleared. + **********************************************************************/ + void geod_polygon_clear(struct geod_polygon* p); + + /** * Add a point to the polygon or polyline. * * @param[in] g a pointer to the geod_geodesic object specifying the @@ -630,6 +754,8 @@ extern "C" { * 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. + * * Example, compute the perimeter and area of the geodesic triangle with * vertices (0°N,0°E), (0°N,90°E), (90°N,0°E). @code{.c} @@ -774,10 +900,7 @@ extern "C" { enum geod_flags { GEOD_NOFLAGS = 0U, /**< No flags */ GEOD_ARCMODE = 1U<<0, /**< Position given in terms of arc distance */ - GEOD_LONG_UNROLL = 1U<<15, /**< Unroll the longitude */ - /**< @cond SKIP */ - GEOD_LONG_NOWRAP = GEOD_LONG_UNROLL /* For backward compatibility only */ - /**< @endcond */ + GEOD_LONG_UNROLL = 1U<<15 /**< Unroll the longitude */ }; #if defined(__cplusplus) diff --git a/src/geodtest.c b/src/geodtest.c new file mode 100644 index 00000000..990a213e --- /dev/null +++ b/src/geodtest.c @@ -0,0 +1,768 @@ +/** + * \file geodtest.c + * \brief Test suite for the geodesic routines in C + * + * Run these tests by configuring with cmake and running "make test". + * + * Copyright (c) Charles Karney (2015-2016) <charles@karney.com> and licensed + * under the MIT/X11 License. For more information, see + * http://geographiclib.sourceforge.net/ + **********************************************************************/ + +/** @cond SKIP */ + +#include "geodesic.h" +#include <stdio.h> +#include <math.h> + +#if defined(_MSC_VER) +// Squelch warnings about assignment within conditional expression +# pragma warning (disable: 4706) +#endif + +double wgs84_a = 6378137, wgs84_f = 1/298.257223563; /* WGS84 */ + +int assertEquals(double x, double y, double d) { + if (fabs(x - y) <= d) + return 0; + printf("assertEquals fails: %.7g != %.7g +/- %.7g\n", x, y, d); + return 1; +} + +const int ncases = 20; +double testcases[20][12] = { + {35.60777, -139.44815, 111.098748429560326, + -11.17491, -69.95921, 129.289270889708762, + 8935244.5604818305, 80.50729714281974, 6273170.2055303837, + 0.16606318447386067, 0.16479116945612937, 12841384694976.432}, + {55.52454, 106.05087, 22.020059880982801, + 77.03196, 197.18234, 109.112041110671519, + 4105086.1713924406, 36.892740690445894, 3828869.3344387607, + 0.80076349608092607, 0.80101006984201008, 61674961290615.615}, + {-21.97856, 142.59065, -32.44456876433189, + 41.84138, 98.56635, -41.84359951440466, + 8394328.894657671, 75.62930491011522, 6161154.5773110616, + 0.24816339233950381, 0.24930251203627892, -6637997720646.717}, + {-66.99028, 112.2363, 173.73491240878403, + -12.70631, 285.90344, 2.512956620913668, + 11150344.2312080241, 100.278634181155759, 6289939.5670446687, + -0.17199490274700385, -0.17722569526345708, -121287239862139.744}, + {-17.42761, 173.34268, -159.033557661192928, + -15.84784, 5.93557, -20.787484651536988, + 16076603.1631180673, 144.640108810286253, 3732902.1583877189, + -0.81273638700070476, -0.81299800519154474, 97825992354058.708}, + {32.84994, 48.28919, 150.492927788121982, + -56.28556, 202.29132, 48.113449399816759, + 16727068.9438164461, 150.565799985466607, 3147838.1910180939, + -0.87334918086923126, -0.86505036767110637, -72445258525585.010}, + {6.96833, 52.74123, 92.581585386317712, + -7.39675, 206.17291, 90.721692165923907, + 17102477.2496958388, 154.147366239113561, 2772035.6169917581, + -0.89991282520302447, -0.89986892177110739, -1311796973197.995}, + {-50.56724, -16.30485, -105.439679907590164, + -33.56571, -94.97412, -47.348547835650331, + 6455670.5118668696, 58.083719495371259, 5409150.7979815838, + 0.53053508035997263, 0.52988722644436602, 41071447902810.047}, + {-58.93002, -8.90775, 140.965397902500679, + -8.91104, 133.13503, 19.255429433416599, + 11756066.0219864627, 105.755691241406877, 6151101.2270708536, + -0.26548622269867183, -0.27068483874510741, -86143460552774.735}, + {-68.82867, -74.28391, 93.774347763114881, + -50.63005, -8.36685, 34.65564085411343, + 3956936.926063544, 35.572254987389284, 3708890.9544062657, + 0.81443963736383502, 0.81420859815358342, -41845309450093.787}, + {-10.62672, -32.0898, -86.426713286747751, + 5.883, -134.31681, -80.473780971034875, + 11470869.3864563009, 103.387395634504061, 6184411.6622659713, + -0.23138683500430237, -0.23155097622286792, 4198803992123.548}, + {-21.76221, 166.90563, 29.319421206936428, + 48.72884, 213.97627, 43.508671946410168, + 9098627.3986554915, 81.963476716121964, 6299240.9166992283, + 0.13965943368590333, 0.14152969707656796, 10024709850277.476}, + {-19.79938, -174.47484, 71.167275780171533, + -11.99349, -154.35109, 65.589099775199228, + 2319004.8601169389, 20.896611684802389, 2267960.8703918325, + 0.93427001867125849, 0.93424887135032789, -3935477535005.785}, + {-11.95887, -116.94513, 92.712619830452549, + 4.57352, 7.16501, 78.64960934409585, + 13834722.5801401374, 124.688684161089762, 5228093.177931598, + -0.56879356755666463, -0.56918731952397221, -9919582785894.853}, + {-87.85331, 85.66836, -65.120313040242748, + 66.48646, 16.09921, -4.888658719272296, + 17286615.3147144645, 155.58592449699137, 2635887.4729110181, + -0.90697975771398578, -0.91095608883042767, 42667211366919.534}, + {1.74708, 128.32011, -101.584843631173858, + -11.16617, 11.87109, -86.325793296437476, + 12942901.1241347408, 116.650512484301857, 5682744.8413270572, + -0.44857868222697644, -0.44824490340007729, 10763055294345.653}, + {-25.72959, -144.90758, -153.647468693117198, + -57.70581, -269.17879, -48.343983158876487, + 9413446.7452453107, 84.664533838404295, 6356176.6898881281, + 0.09492245755254703, 0.09737058264766572, 74515122850712.444}, + {-41.22777, 122.32875, 14.285113402275739, + -7.57291, 130.37946, 10.805303085187369, + 3812686.035106021, 34.34330804743883, 3588703.8812128856, + 0.82605222593217889, 0.82572158200920196, -2456961531057.857}, + {11.01307, 138.25278, 79.43682622782374, + 6.62726, 247.05981, 103.708090215522657, + 11911190.819018408, 107.341669954114577, 6070904.722786735, + -0.29767608923657404, -0.29785143390252321, 17121631423099.696}, + {-29.47124, 95.14681, -163.779130441688382, + -27.46601, -69.15955, -15.909335945554969, + 13487015.8381145492, 121.294026715742277, 5481428.9945736388, + -0.51527225545373252, -0.51556587964721788, 104679964020340.318}}; + +int testinverse() { + double lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12; + double azi1a, azi2a, s12a, a12a, m12a, M12a, M21a, S12a; + struct geod_geodesic g; + int i, result = 0; + geod_init(&g, wgs84_a, wgs84_f); + for (i = 0; i < ncases; ++i) { + lat1 = testcases[i][0]; lon1 = testcases[i][1]; azi1 = testcases[i][2]; + lat2 = testcases[i][3]; lon2 = testcases[i][4]; azi2 = testcases[i][5]; + 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_geninverse(&g, lat1, lon1, lat2, lon2, &s12a, &azi1a, &azi2a, + &m12a, &M12a, &M21a, &S12a); + result += assertEquals(azi1, azi1a, 1e-13); + result += assertEquals(azi2, azi2a, 1e-13); + result += assertEquals(s12, s12a, 1e-8); + result += assertEquals(a12, a12a, 1e-13); + result += assertEquals(m12, m12a, 1e-8); + result += assertEquals(M12, M12a, 1e-15); + result += assertEquals(M21, M21a, 1e-15); + result += assertEquals(S12, S12a, 0.1); + } + return result; +} + +int testdirect() { + double lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12; + double lat2a, lon2a, azi2a, a12a, m12a, M12a, M21a, S12a; + struct geod_geodesic g; + int i, result = 0; + unsigned flags = GEOD_LONG_UNROLL; + geod_init(&g, wgs84_a, wgs84_f); + for (i = 0; i < ncases; ++i) { + lat1 = testcases[i][0]; lon1 = testcases[i][1]; azi1 = testcases[i][2]; + lat2 = testcases[i][3]; lon2 = testcases[i][4]; azi2 = testcases[i][5]; + 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, 0, + &m12a, &M12a, &M21a, &S12a); + result += assertEquals(lat2, lat2a, 1e-13); + result += assertEquals(lon2, lon2a, 1e-13); + result += assertEquals(azi2, azi2a, 1e-13); + result += assertEquals(a12, a12a, 1e-13); + result += assertEquals(m12, m12a, 1e-8); + result += assertEquals(M12, M12a, 1e-15); + result += assertEquals(M21, M21a, 1e-15); + result += assertEquals(S12, S12a, 0.1); + } + return result; +} + +int testarcdirect() { + double lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12; + double lat2a, lon2a, azi2a, s12a, m12a, M12a, M21a, S12a; + struct geod_geodesic g; + int i, result = 0; + unsigned flags = GEOD_ARCMODE | GEOD_LONG_UNROLL; + geod_init(&g, wgs84_a, wgs84_f); + for (i = 0; i < ncases; ++i) { + lat1 = testcases[i][0]; lon1 = testcases[i][1]; azi1 = testcases[i][2]; + lat2 = testcases[i][3]; lon2 = testcases[i][4]; azi2 = testcases[i][5]; + s12 = testcases[i][6]; a12 = testcases[i][7]; m12 = testcases[i][8]; + M12 = testcases[i][9]; M21 = testcases[i][10]; S12 = testcases[i][11]; + geod_gendirect(&g, lat1, lon1, azi1, flags, a12, + &lat2a, &lon2a, &azi2a, &s12a, &m12a, &M12a, &M21a, &S12a); + result += assertEquals(lat2, lat2a, 1e-13); + result += assertEquals(lon2, lon2a, 1e-13); + result += assertEquals(azi2, azi2a, 1e-13); + result += assertEquals(s12, s12a, 1e-8); + result += assertEquals(m12, m12a, 1e-8); + result += assertEquals(M12, M12a, 1e-15); + result += assertEquals(M21, M21a, 1e-15); + result += assertEquals(S12, S12a, 0.1); + } + return result; +} + +int GeodSolve0() { + double azi1, azi2, s12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverse(&g, 40.6, -73.8, 49.01666667, 2.55, &s12, &azi1, &azi2); + result += assertEquals(azi1, 53.47022, 0.5e-5); + result += assertEquals(azi2, 111.59367, 0.5e-5); + result += assertEquals(s12, 5853226, 0.5); + return result; +} + +int GeodSolve1() { + double lat2, lon2, azi2; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_direct(&g, 40.63972222, -73.77888889, 53.5, 5850e3, + &lat2, &lon2, &azi2); + result += assertEquals(lat2, 49.01467, 0.5e-5); + result += assertEquals(lon2, 2.56106, 0.5e-5); + result += assertEquals(azi2, 111.62947, 0.5e-5); + return result; +} + +int GeodSolve2() { + /* Check fix for antipodal prolate bug found 2010-09-04 */ + double azi1, azi2, s12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, 6.4e6, -1/150.0); + geod_inverse(&g, 0.07476, 0, -0.07476, 180, &s12, &azi1, &azi2); + result += assertEquals(azi1, 90.00078, 0.5e-5); + result += assertEquals(azi2, 90.00078, 0.5e-5); + result += assertEquals(s12, 20106193, 0.5); + geod_inverse(&g, 0.1, 0, -0.1, 180, &s12, &azi1, &azi2); + result += assertEquals(azi1, 90.00105, 0.5e-5); + result += assertEquals(azi2, 90.00105, 0.5e-5); + result += assertEquals(s12, 20106193, 0.5); + return result; +} + +int GeodSolve4() { + /* Check fix for short line bug found 2010-05-21 */ + double s12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverse(&g, 36.493349428792, 0, 36.49334942879201, .0000008, + &s12, 0, 0); + result += assertEquals(s12, 0.072, 0.5e-3); + return result; +} + +int GeodSolve5() { + /* Check fix for point2=pole bug found 2010-05-03 */ + double lat2, lon2, azi2; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_direct(&g, 0.01777745589997, 30, 0, 10e6, &lat2, &lon2, &azi2); + result += assertEquals(lat2, 90, 0.5e-5); + if (lon2 < 0) { + result += assertEquals(lon2, -150, 0.5e-5); + result += assertEquals(azi2, -180, 0.5e-5); + } else { + result += assertEquals(lon2, 30, 0.5e-5); + result += assertEquals(azi2, 0, 0.5e-5); + } + return result; +} + +int GeodSolve6() { + /* Check fix for volatile sbet12a bug found 2011-06-25 (gcc 4.4.4 + * x86 -O3). Found again on 2012-03-27 with tdm-mingw32 (g++ 4.6.1). */ + double s12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverse(&g, 88.202499451857, 0, + -88.202499451857, 179.981022032992859592, &s12, 0, 0); + result += assertEquals(s12, 20003898.214, 0.5e-3); + geod_inverse(&g, 89.262080389218, 0, + -89.262080389218, 179.992207982775375662, &s12, 0, 0); + result += assertEquals(s12, 20003925.854, 0.5e-3); + geod_inverse(&g, 89.333123580033, 0, + -89.333123580032997687, 179.99295812360148422, &s12, 0, 0); + result += assertEquals(s12, 20003926.881, 0.5e-3); + return result; +} + +int GeodSolve9() { + /* Check fix for volatile x bug found 2011-06-25 (gcc 4.4.4 x86 -O3) */ + double s12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverse(&g, 56.320923501171, 0, + -56.320923501171, 179.664747671772880215, &s12, 0, 0); + result += assertEquals(s12, 19993558.287, 0.5e-3); + return result; +} + +int GeodSolve10() { + /* Check fix for adjust tol1_ bug found 2011-06-25 (Visual Studio + * 10 rel + debug) */ + double s12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverse(&g, 52.784459512564, 0, + -52.784459512563990912, 179.634407464943777557, &s12, 0, 0); + result += assertEquals(s12, 19991596.095, 0.5e-3); + return result; +} + +int GeodSolve11() { + /* Check fix for bet2 = -bet1 bug found 2011-06-25 (Visual Studio + * 10 rel + debug) */ + double s12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverse(&g, 48.522876735459, 0, + -48.52287673545898293, 179.599720456223079643, &s12, 0, 0); + result += assertEquals(s12, 19989144.774, 0.5e-3); + return result; +} + +int GeodSolve12() { + /* Check fix for inverse geodesics on extreme prolate/oblate + * ellipsoids Reported 2012-08-29 Stefan Guenther + * <stefan.gunther@embl.de>; fixed 2012-10-07 */ + double azi1, azi2, s12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, 89.8, -1.83); + geod_inverse(&g, 0, 0, -10, 160, &s12, &azi1, &azi2); + result += assertEquals(azi1, 120.27, 1e-2); + result += assertEquals(azi2, 105.15, 1e-2); + result += assertEquals(s12, 266.7, 1e-1); + return result; +} + +int GeodSolve14() { + /* Check fix for inverse ignoring lon12 = nan */ + double azi1, azi2, s12, nan = sqrt(-1.0); + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverse(&g, 0, 0, 1, nan, &s12, &azi1, &azi2); + result += azi1 == azi1 ? 1 : 0; + result += azi2 == azi2 ? 1 : 0; + result += s12 == s12 ? 1 : 0; + return result; +} + +int GeodSolve15() { + /* Initial implementation of Math::eatanhe was wrong for e^2 < 0. This + * checks that this is fixed. */ + double S12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, 6.4e6, -1/150.0); + geod_gendirect(&g, 1, 2, 3, 0, 4, + 0, 0, 0, 0, 0, 0, 0, &S12); + result += assertEquals(S12, 23700, 0.5); + return result; +} + +int GeodSolve17() { + /* Check fix for LONG_UNROLL bug found on 2015-05-07 */ + double lat2, lon2, azi2; + struct geod_geodesic g; + struct geod_geodesicline l; + 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, 0, 0, 0, 0, 0); + result += assertEquals(lat2, -39, 1); + result += assertEquals(lon2, -254, 1); + result += assertEquals(azi2, -170, 1); + geod_lineinit(&l, &g, 40, -75, -10, 0); + geod_genposition(&l, flags, 2e7, &lat2, &lon2, &azi2, 0, 0, 0, 0, 0); + result += assertEquals(lat2, -39, 1); + result += assertEquals(lon2, -254, 1); + result += assertEquals(azi2, -170, 1); + geod_direct(&g, 40, -75, -10, 2e7, &lat2, &lon2, &azi2); + result += assertEquals(lat2, -39, 1); + result += assertEquals(lon2, 105, 1); + result += assertEquals(azi2, -170, 1); + geod_position(&l, 2e7, &lat2, &lon2, &azi2); + result += assertEquals(lat2, -39, 1); + result += assertEquals(lon2, 105, 1); + result += assertEquals(azi2, -170, 1); + return result; +} + +int GeodSolve26() { + /* Check 0/0 problem with area calculation on sphere 2015-09-08 */ + double S12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, 6.4e6, 0); + geod_geninverse(&g, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, &S12); + result += assertEquals(S12, 49911046115.0, 0.5); + return result; +} + +int GeodSolve28() { + /* Check for bad placement of assignment of r.a12 with |f| > 0.01 (bug in + * Java implementation fixed on 2015-05-19). */ + double a12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, 6.4e6, 0.1); + a12 = geod_gendirect(&g, 1, 2, 10, 0, 5e6, 0, 0, 0, 0, 0, 0, 0, 0); + result += assertEquals(a12, 48.55570690, 0.5e-8); + return result; +} + +int GeodSolve33() { + /* Check max(-0.0,+0.0) issues 2015-08-22 (triggered by bugs in Octave -- + * sind(-0.0) = +0.0 -- and in some version of Visual Studio -- + * fmod(-0.0, 360.0) = +0.0. */ + double azi1, azi2, s12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverse(&g, 0, 0, 0, 179, &s12, &azi1, &azi2); + result += assertEquals(azi1, 90.00000, 0.5e-5); + result += assertEquals(azi2, 90.00000, 0.5e-5); + result += assertEquals(s12, 19926189, 0.5); + geod_inverse(&g, 0, 0, 0, 179.5, &s12, &azi1, &azi2); + result += assertEquals(azi1, 55.96650, 0.5e-5); + result += assertEquals(azi2, 124.03350, 0.5e-5); + result += assertEquals(s12, 19980862, 0.5); + geod_inverse(&g, 0, 0, 0, 180, &s12, &azi1, &azi2); + result += assertEquals(azi1, 0.00000, 0.5e-5); + result += assertEquals(azi2, -180.00000, 0.5e-5); + result += assertEquals(s12, 20003931, 0.5); + geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2); + result += assertEquals(azi1, 0.00000, 0.5e-5); + result += assertEquals(azi2, -180.00000, 0.5e-5); + result += assertEquals(s12, 19893357, 0.5); + geod_init(&g, 6.4e6, 0); + geod_inverse(&g, 0, 0, 0, 179, &s12, &azi1, &azi2); + result += assertEquals(azi1, 90.00000, 0.5e-5); + result += assertEquals(azi2, 90.00000, 0.5e-5); + result += assertEquals(s12, 19994492, 0.5); + geod_inverse(&g, 0, 0, 0, 180, &s12, &azi1, &azi2); + result += assertEquals(azi1, 0.00000, 0.5e-5); + result += assertEquals(azi2, -180.00000, 0.5e-5); + result += assertEquals(s12, 20106193, 0.5); + geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2); + result += assertEquals(azi1, 0.00000, 0.5e-5); + result += assertEquals(azi2, -180.00000, 0.5e-5); + result += assertEquals(s12, 19994492, 0.5); + geod_init(&g, 6.4e6, -1/300.0); + geod_inverse(&g, 0, 0, 0, 179, &s12, &azi1, &azi2); + result += assertEquals(azi1, 90.00000, 0.5e-5); + result += assertEquals(azi2, 90.00000, 0.5e-5); + result += assertEquals(s12, 19994492, 0.5); + geod_inverse(&g, 0, 0, 0, 180, &s12, &azi1, &azi2); + result += assertEquals(azi1, 90.00000, 0.5e-5); + result += assertEquals(azi2, 90.00000, 0.5e-5); + result += assertEquals(s12, 20106193, 0.5); + geod_inverse(&g, 0, 0, 0.5, 180, &s12, &azi1, &azi2); + result += assertEquals(azi1, 33.02493, 0.5e-5); + result += assertEquals(azi2, 146.97364, 0.5e-5); + result += assertEquals(s12, 20082617, 0.5); + geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2); + result += assertEquals(azi1, 0.00000, 0.5e-5); + result += assertEquals(azi2, -180.00000, 0.5e-5); + result += assertEquals(s12, 20027270, 0.5); + + return result; +} + +int GeodSolve55() { + /* Check fix for nan + point on equator or pole not returning all nans in + * Geodesic::Inverse, found 2015-09-23. */ + double azi1, azi2, s12, nan = sqrt(-1.0); + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverse(&g, nan, 0, 0, 90, &s12, &azi1, &azi2); + result += azi1 == azi1 ? 1 : 0; + result += azi2 == azi2 ? 1 : 0; + result += s12 == s12 ? 1 : 0; + geod_inverse(&g, nan, 0, 90, 9, &s12, &azi1, &azi2); + result += azi1 == azi1 ? 1 : 0; + result += azi2 == azi2 ? 1 : 0; + result += s12 == s12 ? 1 : 0; + return result; +} + +int GeodSolve59() { + /* Check for points close with longitudes close to 180 deg apart. */ + double azi1, azi2, s12; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverse(&g, 5, 0.00000000000001, 10, 180, &s12, &azi1, &azi2); + result += assertEquals(azi1, 0.000000000000035, 1.5e-14); + result += assertEquals(azi2, 179.99999999999996, 1.5e-14); + result += assertEquals(s12, 18345191.174332713, 2.5e-9); + return result; +} + +int GeodSolve61() { + /* Make sure small negative azimuths are west-going */ + double lat2, lon2, azi2; + struct geod_geodesic g; + struct geod_geodesicline l; + int result = 0; + unsigned flags = GEOD_LONG_UNROLL; + geod_init(&g, wgs84_a, wgs84_f); + geod_gendirect(&g, 45, 0, -0.000000000000000003, flags, 1e7, + &lat2, &lon2, &azi2, 0, 0, 0, 0, 0); + result += assertEquals(lat2, 45.30632, 0.5e-5); + result += assertEquals(lon2, -180, 0.5e-5); + result += assertEquals(azi2, -180, 0.5e-5); + geod_inverseline(&l, &g, 45, 0, 80, -0.000000000000000003, 0); + geod_genposition(&l, flags, 1e7, &lat2, &lon2, &azi2, 0, 0, 0, 0, 0); + result += assertEquals(lat2, 45.30632, 0.5e-5); + result += assertEquals(lon2, -180, 0.5e-5); + result += assertEquals(azi2, -180, 0.5e-5); + return result; +} + +int GeodSolve65() { + /* Check for bug in east-going check in GeodesicLine (needed to check for + * sign of 0) and sign error in area calculation due to a bogus override of + * the code for alp12. Found/fixed on 2015-12-19. */ + double lat2, lon2, azi2, s12, a12, m12, M12, M21, S12; + struct geod_geodesic g; + struct geod_geodesicline l; + int result = 0; + unsigned flags = GEOD_LONG_UNROLL, caps = GEOD_ALL; + geod_init(&g, wgs84_a, wgs84_f); + geod_inverseline(&l, &g, 30, -0.000000000000000001, -31, 180, caps); + a12 = geod_genposition(&l, flags, 1e7, + &lat2, &lon2, &azi2, &s12, &m12, &M12, &M21, &S12); + result += assertEquals(lat2, -60.23169, 0.5e-5); + result += assertEquals(lon2, -0.00000, 0.5e-5); + result += assertEquals(azi2, -180.00000, 0.5e-5); + result += assertEquals(s12, 10000000, 0.5); + result += assertEquals(a12, 90.06544, 0.5e-5); + result += assertEquals(m12, 6363636, 0.5); + result += assertEquals(M12, -0.0012834, 0.5e-7); + result += assertEquals(M21, 0.0013749, 0.5e-7); + result += assertEquals(S12, 0, 0.5); + a12 = geod_genposition(&l, flags, 2e7, + &lat2, &lon2, &azi2, &s12, &m12, &M12, &M21, &S12); + result += assertEquals(lat2, -30.03547, 0.5e-5); + result += assertEquals(lon2, -180.00000, 0.5e-5); + result += assertEquals(azi2, -0.00000, 0.5e-5); + result += assertEquals(s12, 20000000, 0.5); + result += assertEquals(a12, 179.96459, 0.5e-5); + result += assertEquals(m12, 54342, 0.5); + result += assertEquals(M12, -1.0045592, 0.5e-7); + result += assertEquals(M21, -0.9954339, 0.5e-7); + result += assertEquals(S12, 127516405431022.0, 0.5); + return result; +} + +int GeodSolve67() { + /* Check for InverseLine if line is slightly west of S and that s13 is + correctly set. */ + double lat2, lon2, azi2; + struct geod_geodesic g; + struct geod_geodesicline l; + int result = 0; + 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, 0, 0, 0, 0, 0); + result += assertEquals(lat2, 4.96445, 0.5e-5); + result += assertEquals(lon2, -180.00000, 0.5e-5); + result += assertEquals(azi2, -0.00000, 0.5e-5); + geod_genposition(&l, flags, 0.5 * l.s13, &lat2, &lon2, &azi2, 0, 0, 0, 0, 0); + result += assertEquals(lat2, -87.52461, 0.5e-5); + result += assertEquals(lon2, -0.00000, 0.5e-5); + result += assertEquals(azi2, -180.00000, 0.5e-5); + return result; +} + +int GeodSolve71() { + /* Check that DirectLine sets s13. */ + double lat2, lon2, azi2; + struct geod_geodesic g; + struct geod_geodesicline l; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_directline(&l, &g, 1, 2, 45, 1e7, 0); + geod_position(&l, 0.5 * l.s13, &lat2, &lon2, &azi2); + result += assertEquals(lat2, 30.92625, 0.5e-5); + result += assertEquals(lon2, 37.54640, 0.5e-5); + result += assertEquals(azi2, 55.43104, 0.5e-5); + return result; +} + +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. */ + double lat2, lon2, azi2; + struct geod_geodesic g; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_direct(&g, 90, 10, 180, -1e6, + &lat2, &lon2, &azi2); + result += assertEquals(lat2, 81.04623, 0.5e-5); + result += assertEquals(lon2, -170, 0.5e-5); + result += assertEquals(azi2, 0, 0.5e-5); + return result; +} + +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); + for (i = 0; i < N; ++i) + geod_polygon_addpoint(g, &p, points[i][0], points[i][1]); + geod_polygon_compute(g, &p, 0, 1, area, perimeter); +} + +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); + for (i = 0; i < N; ++i) + geod_polygon_addpoint(g, &p, points[i][0], points[i][1]); + geod_polygon_compute(g, &p, 0, 1, 0, perimeter); +} + +int Planimeter0() { + /* Check fix for pole-encircling bug found 2011-03-16 */ + double pa[4][2] = {{89, 0}, {89, 90}, {89, 180}, {89, 270}}; + double pb[4][2] = {{-89, 0}, {-89, 90}, {-89, 180}, {-89, 270}}; + double pc[4][2] = {{0, -1}, {-1, 0}, {0, 1}, {1, 0}}; + double pd[3][2] = {{90, 0}, {0, 0}, {0, 90}}; + struct geod_geodesic g; + double perimeter, area; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + + planimeter(&g, pa, 4, &perimeter, &area); + result += assertEquals(perimeter, 631819.8745, 1e-4); + result += assertEquals(area, 24952305678.0, 1); + + planimeter(&g, pb, 4, &perimeter, &area); + result += assertEquals(perimeter, 631819.8745, 1e-4); + result += assertEquals(area, -24952305678.0, 1); + + planimeter(&g, pc, 4, &perimeter, &area); + result += assertEquals(perimeter, 627598.2731, 1e-4); + result += assertEquals(area, 24619419146.0, 1); + + planimeter(&g, pd, 3, &perimeter, &area); + result += assertEquals(perimeter, 30022685, 1); + result += assertEquals(area, 63758202715511.0, 1); + + polylength(&g, pd, 3, &perimeter); + result += assertEquals(perimeter, 20020719, 1); + + return result; +} + +int Planimeter5() { + /* Check fix for Planimeter pole crossing bug found 2011-06-24 */ + double points[3][2] = {{89, 0.1}, {89, 90.1}, {89, -179.9}}; + struct geod_geodesic g; + double perimeter, area; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + planimeter(&g, points, 3, &perimeter, &area); + result += assertEquals(perimeter, 539297, 1); + result += assertEquals(area, 12476152838.5, 1); + return result; +} + +int Planimeter6() { + /* Check fix for Planimeter lon12 rounding bug found 2012-12-03 */ + double pa[3][2] = {{9, -0.00000000000001}, {9, 180}, {9, 0}}; + double pb[3][2] = {{9, 0.00000000000001}, {9, 0}, {9, 180}}; + double pc[3][2] = {{9, 0.00000000000001}, {9, 180}, {9, 0}}; + double pd[3][2] = {{9, -0.00000000000001}, {9, 0}, {9, 180}}; + struct geod_geodesic g; + double perimeter, area; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + + planimeter(&g, pa, 3, &perimeter, &area); + result += assertEquals(perimeter, 36026861, 1); + result += assertEquals(area, 0, 1); + planimeter(&g, pb, 3, &perimeter, &area); + result += assertEquals(perimeter, 36026861, 1); + result += assertEquals(area, 0, 1); + planimeter(&g, pc, 3, &perimeter, &area); + result += assertEquals(perimeter, 36026861, 1); + result += assertEquals(area, 0, 1); + planimeter(&g, pd, 3, &perimeter, &area); + result += assertEquals(perimeter, 36026861, 1); + result += assertEquals(area, 0, 1); + return result; +} + +int Planimeter12() { + /* Area of arctic circle (not really -- adjunct to rhumb-area test) */ + double points[2][2] = {{66.562222222, 0}, {66.562222222, 180}}; + struct geod_geodesic g; + double perimeter, area; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + planimeter(&g, points, 2, &perimeter, &area); + result += assertEquals(perimeter, 10465729, 1); + result += assertEquals(area, 0, 1); + return result; +} + +int Planimeter13() { + /* Check encircling pole twice */ + double points[6][2] = {{89,-360}, {89,-240}, {89,-120}, + {89,0}, {89,120}, {89,240}}; + struct geod_geodesic g; + double perimeter, area; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + planimeter(&g, points, 6, &perimeter, &area); + result += assertEquals(perimeter, 1160741, 1); + result += assertEquals(area, 32415230256.0, 1); + return result; +} + +int main() { + int n = 0, i; + if ((i = testinverse())) {++n; printf("testinverse fail: %d\n", i);} + if ((i = testdirect())) {++n; printf("testdirect fail: %d\n", i);} + if ((i = testarcdirect())) {++n; printf("testarcdirect fail: %d\n", i);} + if ((i = GeodSolve0())) {++n; printf("GeodSolve0 fail: %d\n", i);} + if ((i = GeodSolve1())) {++n; printf("GeodSolve1 fail: %d\n", i);} + if ((i = GeodSolve2())) {++n; printf("GeodSolve2 fail: %d\n", i);} + if ((i = GeodSolve4())) {++n; printf("GeodSolve4 fail: %d\n", i);} + if ((i = GeodSolve5())) {++n; printf("GeodSolve5 fail: %d\n", i);} + if ((i = GeodSolve6())) {++n; printf("GeodSolve6 fail: %d\n", i);} + if ((i = GeodSolve9())) {++n; printf("GeodSolve9 fail: %d\n", i);} + if ((i = GeodSolve10())) {++n; printf("GeodSolve10 fail: %d\n", i);} + if ((i = GeodSolve11())) {++n; printf("GeodSolve11 fail: %d\n", i);} + if ((i = GeodSolve12())) {++n; printf("GeodSolve12 fail: %d\n", i);} + if ((i = GeodSolve14())) {++n; printf("GeodSolve14 fail: %d\n", i);} + if ((i = GeodSolve15())) {++n; printf("GeodSolve15 fail: %d\n", i);} + if ((i = GeodSolve17())) {++n; printf("GeodSolve17 fail: %d\n", i);} + if ((i = GeodSolve26())) {++n; printf("GeodSolve26 fail: %d\n", i);} + if ((i = GeodSolve28())) {++n; printf("GeodSolve28 fail: %d\n", i);} + if ((i = GeodSolve33())) {++n; printf("GeodSolve33 fail: %d\n", i);} + if ((i = GeodSolve55())) {++n; printf("GeodSolve55 fail: %d\n", i);} + if ((i = GeodSolve59())) {++n; printf("GeodSolve59 fail: %d\n", i);} + if ((i = GeodSolve61())) {++n; printf("GeodSolve61 fail: %d\n", i);} + if ((i = GeodSolve65())) {++n; printf("GeodSolve65 fail: %d\n", i);} + if ((i = GeodSolve67())) {++n; printf("GeodSolve67 fail: %d\n", i);} + if ((i = GeodSolve71())) {++n; printf("GeodSolve71 fail: %d\n", i);} + if ((i = GeodSolve73())) {++n; printf("GeodSolve73 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);} + if ((i = Planimeter12())) {++n; printf("Planimeter12 fail: %d\n", i);} + if ((i = Planimeter13())) {++n; printf("Planimeter13 fail: %d\n", i);} + return n; +} + +/** @endcond */ diff --git a/src/pj_release.c b/src/pj_release.c index c3736b44..1f5a201e 100644 --- a/src/pj_release.c +++ b/src/pj_release.c @@ -2,7 +2,7 @@ #include <projects.h> -char const pj_release[]="Rel. 4.9.2, 08 September 2015"; +char const pj_release[]="Rel. 4.9.3, dd Month yyyy"; const char *pj_get_release() diff --git a/src/proj_api.h b/src/proj_api.h index e381815c..db833b11 100644 --- a/src/proj_api.h +++ b/src/proj_api.h @@ -38,7 +38,7 @@ extern "C" { #endif /* Try to update this every version! */ -#define PJ_VERSION 492 +#define PJ_VERSION 493 /* pj_init() and similar functions can be used with a non-C locale */ /* Can be detected too at runtime if the symbol pj_atof exists */ |