diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/geodesic.c | 86 | ||||
| -rw-r--r-- | src/geodesic.h | 41 |
2 files changed, 85 insertions, 42 deletions
diff --git a/src/geodesic.c b/src/geodesic.c index c5d481ee..957f2144 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -195,6 +195,8 @@ static real InverseStart(const struct geod_geodesic* g, real* psalp1, real* pcalp1, /* Only updated if return val >= 0 */ real* psalp2, real* pcalp2, + /* Only updated for short lines */ + real* pdnm, /* Scratch areas of the right size */ real C1a[], real C2a[]); static real Lambda12(const struct geod_geodesic* g, @@ -218,6 +220,8 @@ static void C1pf(real eps, real c[]); static real A2m1f(real eps); static void C2f(real eps, real c[]); static int transit(real lon1, real lon2); +static void accini(real s[]); +static void accadd(real s[], real y); /** @endcond */ @@ -234,8 +238,18 @@ void geod_init(struct geod_geodesic* g, real a, real f) { (g->e2 == 0 ? 1 : (g->e2 > 0 ? atanhx(sqrt(g->e2)) : atan(sqrt(-g->e2))) / sqrt(fabs(g->e2))))/2; /* authalic radius squared */ - /* The sig12 threshold for "really short" */ - g->etol2 = 0.01 * tol2 / maxx((real)(0.1), sqrt(fabs(g->e2))); + /* The sig12 threshold for "really short". Using the auxiliary sphere + * solution with dnm computed at (bet1 + bet2) / 2, the relative error in the + * azimuth consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2. (Error + * measured for 1/100 < b/a < 100 and abs(f) >= 1/1000. For a given f and + * sig12, the max error occurs for lines near the pole. If the old rule for + * computing dnm = (dn1 + dn2)/2 is used, then the error increases by a + * factor of 2.) Setting this equal to epsilon gives sig12 = etol2. Here + * 0.1 is a safety factor (error decreased by 100) and max(0.001, abs(f)) + * stops etol2 getting too large in the nearly spherical case. */ + g->etol2 = 0.1 * tol2 / + sqrt( maxx((real)(0.001), fabs(g->f)) * minx((real)(1), 1 - g->f/2) / 2 ); + A3coeff(g); C3coeff(g); C4coeff(g); @@ -716,14 +730,14 @@ real geod_geninverse(const struct geod_geodesic* g, * meridian and geodesic is neither meridional or equatorial. */ /* Figure a starting point for Newton's method */ + real dnm = 0; sig12 = InverseStart(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12, - &salp1, &calp1, &salp2, &calp2, + &salp1, &calp1, &salp2, &calp2, &dnm, C1a, C2a); if (sig12 >= 0) { - /* Short lines (InverseStart sets salp2, calp2) */ - real dnm = (dn1 + dn2) / 2; + /* Short lines (InverseStart sets salp2, calp2, dnm) */ s12x = sig12 * g->b * dnm; m12x = sq(dnm) * g->b * sin(sig12 / dnm); if (outmask & GEOD_GEODESICSCALE) @@ -1046,9 +1060,11 @@ real InverseStart(const struct geod_geodesic* g, real* psalp1, real* pcalp1, /* Only updated if return val >= 0 */ real* psalp2, real* pcalp2, + /* Only updated for short lines */ + real* pdnm, /* Scratch areas of the right size */ real C1a[], real C2a[]) { - real salp1 = 0, calp1 = 0, salp2 = 0, calp2 = 0; + real salp1 = 0, calp1 = 0, salp2 = 0, calp2 = 0, dnm = 0; /* Return a starting point for Newton's method in salp1 and calp1 (function * value is -1). If Newton's method doesn't need to be used, return also @@ -1076,11 +1092,17 @@ real InverseStart(const struct geod_geodesic* g, real sbet12a = sbet2 * cbet1 + cbet2 * sbet1; #endif boolx shortline = cbet12 >= 0 && sbet12 < (real)(0.5) && - lam12 <= pi / 6; - real - omg12 = !shortline ? lam12 : lam12 / (g->f1 * (dn1 + dn2) / 2), - somg12 = sin(omg12), comg12 = cos(omg12); - real ssig12, csig12; + cbet2 * lam12 < (real)(0.5); + real omg12 = lam12, somg12, comg12, ssig12, csig12; + if (shortline) { + real sbetm2 = sq(sbet1 + sbet2); + /* 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; + } + somg12 = sin(omg12); comg12 = cos(omg12); salp1 = cbet2 * somg12; calp1 = comg12 >= 0 ? @@ -1093,7 +1115,8 @@ real InverseStart(const struct geod_geodesic* g, if (shortline && ssig12 < g->etol2) { /* really short lines */ salp2 = cbet1 * somg12; - calp2 = sbet12 - cbet1 * sbet2 * sq(somg12) / (1 + comg12); + calp2 = sbet12 - cbet1 * sbet2 * + (comg12 >= 0 ? sq(somg12) / (1 + comg12) : 1 - comg12); SinCosNorm(&salp2, &calp2); /* Set return value */ sig12 = atan2(ssig12, csig12); @@ -1200,6 +1223,8 @@ real InverseStart(const struct geod_geodesic* g, *psalp1 = salp1; *pcalp1 = calp1; + if (shortline) + *pdnm = dnm; if (sig12 >= 0) { *psalp2 = salp2; *pcalp2 = calp2; @@ -1492,28 +1517,45 @@ int transit(real lon1, real lon2) { (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0); } +void accini(real s[]) { + /* Initialize an accumulator; this is an array with two elements. */ + s[0] = s[1] = 0; +} + +void accadd(real s[], real y) { + /* Add y to an accumulator. */ + real u, z = sumx(y, s[1], &u); + s[0] = sumx(z, s[0], &s[1]); + if (s[0] == 0) + s[0] = u; + else + s[1] = s[1] + u; +} + /** @endcond */ void geod_polygonarea(const struct geod_geodesic* g, real lats[], real lons[], int n, real* pA, real* pP) { int i, crossings = 0; - real area0 = 4 * pi * g->c2, A = 0, P = 0; + real area0 = 4 * pi * g->c2, A[2], P[2]; + accini(A); accini(P); for (i = 0; i < n; ++i) { real s12, S12; geod_geninverse(g, lats[i], lons[i], lats[(i + 1) % n], lons[(i + 1) % n], &s12, 0, 0, 0, 0, 0, &S12); - P += s12; - A -= S12; /* The minus sign is due to the counter-clockwise convention */ + accadd(P, s12); + /* The minus sign is due to the counter-clockwise convention */ + accadd(A, -S12); crossings += transit(lons[i], lons[(i + 1) % n]); } if (crossings & 1) - A += (A < 0 ? 1 : -1) * area0/2; + accadd(A, (A[0] < 0 ? 1 : -1) * area0/2); /* Put area in (-area0/2, area0/2] */ - if (A > area0/2) - A -= area0; - else if (A <= -area0/2) - A += area0; - if (pA) *pA = A; - if (pP) *pP = P; + if (A[0] > area0/2) + accadd(A, -area0); + else if (A[0] <= -area0/2) + accadd(A, +area0); + if (pA) *pA = A[0]; + if (pP) *pP = P[0]; } diff --git a/src/geodesic.h b/src/geodesic.h index 2c49d0dc..87d087a5 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -116,7 +116,7 @@ * http://geographiclib.sourceforge.net/ * * This library was distributed with - * <a href="../index.html">GeographicLib</a> 1.30. + * <a href="../index.html">GeographicLib</a> 1.31. **********************************************************************/ #if !defined(GEODESIC_H) @@ -127,7 +127,8 @@ extern "C" { #endif /** - * The struct containing information about the ellipsoid. + * The struct containing information about the ellipsoid. This must be + * initialized by geod_init() before use. **********************************************************************/ struct geod_geodesic { double a; /**< the equatorial radius */ @@ -139,7 +140,8 @@ extern "C" { }; /** - * The struct containing information about a single geodesic. + * The struct containing information about a single geodesic. This must be + * initialized by geod_lineinit() before use. **********************************************************************/ struct geod_geodesicline { double lat1; /**< the starting latitude */ @@ -182,17 +184,17 @@ extern "C" { * should be in the range [−90°, 90°]; \e azi1 should be in the * range [−540°, 540°). * - * The geod_mask values are + * 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 + * 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 + * 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 + * 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. @@ -226,11 +228,11 @@ extern "C" { * need some quantities computed. * * If either point is at a pole, the azimuth is defined by keeping the - * longitude fixed and writing \e lat = ±(90° − ε) - * and taking the limit ε → 0+. An arc length greater that - * 180° signifies a geodesic which is not a shortest path. (For a - * prolate ellipsoid, an additional condition is necessary for a shortest - * path: the longitudinal extent must not exceed of 180°.) + * longitude fixed, writing \e lat = ±(90° − ε), and + * taking the limit ε → 0+. An arc length greater that 180° + * signifies a geodesic which is not a shortest path. (For a prolate + * ellipsoid, an additional condition is necessary for a shortest path: the + * longitudinal extent must not exceed of 180°.) * * Example, determine the point 10000 km NE of JFK: @code @@ -266,8 +268,8 @@ extern "C" { * not need some quantities computed. * * If either point is at a pole, the azimuth is defined by keeping the - * longitude fixed and writing \e lat = ±(90° − ε) - * and taking the limit ε → 0+. + * longitude fixed, writing \e lat = ±(90° − ε), and + * taking the limit ε → 0+. * * The solution to the inverse problem is found using Newton's method. If * this fails to converge (this is very unlikely in geodetic applications @@ -436,7 +438,6 @@ 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). * * \e l must have been initialized with a call to geod_lineinit() with \e @@ -514,7 +515,7 @@ extern "C" { * mask values for geod_geodesicline capabilities. **********************************************************************/ enum geod_mask { - GEOD_NONE = 0U, + GEOD_NONE = 0U, /**< Calculate nothing */ GEOD_LATITUDE = 1U<<7 | 0U, /**< Calculate latitude */ GEOD_LONGITUDE = 1U<<8 | 1U<<3, /**< Calculate longitude */ GEOD_AZIMUTH = 1U<<9 | 0U, /**< Calculate azimuth */ |