From 9a55b4f828713e5d109b30afb01d4120e3b25c15 Mon Sep 17 00:00:00 2001 From: Alexandre Barreto Date: Thu, 5 Jan 2023 13:44:16 -0500 Subject: [PATCH] install.rst: improve instructions regarding proj-data (fixes #3539) --- docs/source/install.rst | 4 +-- src/geodesic.c | 62 ++++++++++++++++++++++------------------ src/geodesic.h | 2 +- src/init.cpp | 2 +- src/tests/geodsigntest.c | 19 ++++++++++++ 5 files changed, 57 insertions(+), 32 deletions(-) diff --git a/docs/source/install.rst b/docs/source/install.rst index edd1d32dfaa2..f23e8e98fcaa 100644 --- a/docs/source/install.rst +++ b/docs/source/install.rst @@ -206,8 +206,8 @@ can be modified to suit the users needs. See :ref:`projsync` for more options. As an alternative on systems where network access is disabled, the :ref:`proj-data ` - package can be downloaded and added to the :envvar:`PROJ_DATA` directory - (called ``PROJ_LIB`` before PROJ 9.1) + package can be downloaded and its content decompressed into one of the + directories where PROJ looks for :ref:`resources ` Starting with PROJ 9.2, a ``uninstall`` target is available to remove files installed by the ``install`` target:: diff --git a/src/geodesic.c b/src/geodesic.c index d86e508fc183..6b5a4792eb38 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -77,7 +77,7 @@ static void Init(void) { tol1 = 200 * tol0; tol2 = sqrt(tol0); /* Check on bisection interval */ - tolb = tol0 * tol2; + tolb = tol0; xthresh = 1000 * tol2; degree = pi/hd; NaN = nan("0"); @@ -111,7 +111,7 @@ static double sumx(double u, double v, double* t) { return s; } -static double polyval(int N, const double p[], double x) { +static double polyvalx(int N, const double p[], double x) { double y = N < 0 ? 0 : *p++; while (--N >= 0) y = y * x + *p++; return y; @@ -879,7 +879,7 @@ static double geod_geninverse_int(const struct geod_geodesic* g, double salp1a = tiny, calp1a = 1, salp1b = tiny, calp1b = -1; boolx tripn = FALSE; boolx tripb = FALSE; - for (; numit < maxit2; ++numit) { + for (;; ++numit) { /* the WGS84 test set: mean = 1.47, sd = 1.25, max = 16 * WGS84 and random input: mean = 2.85, sd = 0.60 */ double dv = 0, @@ -887,8 +887,12 @@ static double geod_geninverse_int(const struct geod_geodesic* g, slam12, clam12, &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2, &eps, &domg12, numit < maxit1, &dv, Ca); - /* Reversed test to allow escape with NaNs */ - if (tripb || !(fabs(v) >= (tripn ? 8 : 1) * tol0)) break; + if (tripb || + /* Reversed test to allow escape with NaNs */ + !(fabs(v) >= (tripn ? 8 : 1) * tol0) || + /* Enough bisections to get accurate result */ + numit == maxit2) + break; /* Update bracketing values */ if (v > 0 && (numit > maxit1 || calp1/salp1 > calp1b/salp1b)) { salp1b = salp1; calp1b = calp1; } @@ -897,18 +901,20 @@ static double geod_geninverse_int(const struct geod_geodesic* g, if (numit < maxit1 && dv > 0) { double dalp1 = -v/dv; - double - sdalp1 = sin(dalp1), cdalp1 = cos(dalp1), - nsalp1 = salp1 * cdalp1 + calp1 * sdalp1; - if (nsalp1 > 0 && fabs(dalp1) < pi) { - calp1 = calp1 * cdalp1 - salp1 * sdalp1; - salp1 = nsalp1; - norm2(&salp1, &calp1); - /* In some regimes we don't get quadratic convergence because - * slope -> 0. So use convergence conditions based on epsilon - * instead of sqrt(epsilon). */ - tripn = fabs(v) <= 16 * tol0; - continue; + if (fabs(dalp1) < pi) { + double + sdalp1 = sin(dalp1), cdalp1 = cos(dalp1), + nsalp1 = salp1 * cdalp1 + calp1 * sdalp1; + if (nsalp1 > 0) { + calp1 = calp1 * cdalp1 - salp1 * sdalp1; + salp1 = nsalp1; + norm2(&salp1, &calp1); + /* In some regimes we don't get quadratic convergence because + * slope -> 0. So use convergence conditions based on epsilon + * instead of sqrt(epsilon). */ + tripn = fabs(v) <= 16 * tol0; + continue; + } } } /* Either dv was not positive or updated value was outside legal @@ -1480,7 +1486,7 @@ double Lambda12(const struct geod_geodesic* g, double A3f(const struct geod_geodesic* g, double eps) { /* Evaluate A3 */ - return polyval(nA3 - 1, g->A3x, eps); + return polyvalx(nA3 - 1, g->A3x, eps); } void C3f(const struct geod_geodesic* g, double eps, double c[]) { @@ -1491,7 +1497,7 @@ void C3f(const struct geod_geodesic* g, double eps, double c[]) { for (l = 1; l < nC3; ++l) { /* l is index of C3[l] */ int m = nC3 - l - 1; /* order of polynomial in eps */ mult *= eps; - c[l] = mult * polyval(m, g->C3x + o, eps); + c[l] = mult * polyvalx(m, g->C3x + o, eps); o += m + 1; } } @@ -1503,7 +1509,7 @@ void C4f(const struct geod_geodesic* g, double eps, double c[]) { int o = 0, l; for (l = 0; l < nC4; ++l) { /* l is index of C4[l] */ int m = nC4 - l - 1; /* order of polynomial in eps */ - c[l] = mult * polyval(m, g->C4x + o, eps); + c[l] = mult * polyvalx(m, g->C4x + o, eps); o += m + 1; mult *= eps; } @@ -1516,7 +1522,7 @@ double A1m1f(double eps) { 1, 4, 64, 0, 256, }; int m = nA1/2; - double t = polyval(m, coeff, sq(eps)) / coeff[m + 1]; + double t = polyvalx(m, coeff, sq(eps)) / coeff[m + 1]; return (t + eps) / (1 - eps); } @@ -1542,7 +1548,7 @@ void C1f(double eps, double c[]) { int o = 0, l; for (l = 1; l <= nC1; ++l) { /* l is index of C1p[l] */ int m = (nC1 - l) / 2; /* order of polynomial in eps^2 */ - c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1]; + c[l] = d * polyvalx(m, coeff + o, eps2) / coeff[o + m + 1]; o += m + 2; d *= eps; } @@ -1570,7 +1576,7 @@ void C1pf(double eps, double c[]) { int o = 0, l; for (l = 1; l <= nC1p; ++l) { /* l is index of C1p[l] */ int m = (nC1p - l) / 2; /* order of polynomial in eps^2 */ - c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1]; + c[l] = d * polyvalx(m, coeff + o, eps2) / coeff[o + m + 1]; o += m + 2; d *= eps; } @@ -1583,7 +1589,7 @@ double A2m1f(double eps) { -11, -28, -192, 0, 256, }; int m = nA2/2; - double t = polyval(m, coeff, sq(eps)) / coeff[m + 1]; + double t = polyvalx(m, coeff, sq(eps)) / coeff[m + 1]; return (t - eps) / (1 + eps); } @@ -1609,7 +1615,7 @@ void C2f(double eps, double c[]) { int o = 0, l; for (l = 1; l <= nC2; ++l) { /* l is index of C2[l] */ int m = (nC2 - l) / 2; /* order of polynomial in eps^2 */ - c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1]; + c[l] = d * polyvalx(m, coeff + o, eps2) / coeff[o + m + 1]; o += m + 2; d *= eps; } @@ -1634,7 +1640,7 @@ void A3coeff(struct geod_geodesic* g) { int o = 0, k = 0, j; for (j = nA3 - 1; j >= 0; --j) { /* coeff of eps^j */ int m = nA3 - j - 1 < j ? nA3 - j - 1 : j; /* order of polynomial in n */ - g->A3x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1]; + g->A3x[k++] = polyvalx(m, coeff + o, g->n) / coeff[o + m + 1]; o += m + 2; } } @@ -1677,7 +1683,7 @@ void C3coeff(struct geod_geodesic* g) { for (l = 1; l < nC3; ++l) { /* l is index of C3[l] */ for (j = nC3 - 1; j >= l; --j) { /* coeff of eps^j */ int m = nC3 - j - 1 < j ? nC3 - j - 1 : j; /* order of polynomial in n */ - g->C3x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1]; + g->C3x[k++] = polyvalx(m, coeff + o, g->n) / coeff[o + m + 1]; o += m + 2; } } @@ -1733,7 +1739,7 @@ void C4coeff(struct geod_geodesic* g) { for (l = 0; l < nC4; ++l) { /* l is index of C4[l] */ for (j = nC4 - 1; j >= l; --j) { /* coeff of eps^j */ int m = nC4 - j - 1; /* order of polynomial in n */ - g->C4x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1]; + g->C4x[k++] = polyvalx(m, coeff + o, g->n) / coeff[o + m + 1]; o += m + 2; } } diff --git a/src/geodesic.h b/src/geodesic.h index 8ecb771a4501..db4bf0dafa42 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -31,7 +31,7 @@ * The minor version of the geodesic library. (This tracks the version of * GeographicLib.) **********************************************************************/ -#define GEODESIC_VERSION_MINOR 0 +#define GEODESIC_VERSION_MINOR 1 /** * The patch level of the geodesic library. (This tracks the version of * GeographicLib.) diff --git a/src/init.cpp b/src/init.cpp index fd4a7b2f97b6..4d73f0124dec 100644 --- a/src/init.cpp +++ b/src/init.cpp @@ -772,7 +772,7 @@ pj_init_ctx_with_allow_init_epsg(PJ_CONTEXT *ctx, int argc, char **argv, int all PIN->geod = static_cast(calloc (1, sizeof (struct geod_geodesic))); if (nullptr==PIN->geod) return pj_default_destructor (PIN, PROJ_ERR_OTHER /*ENOMEM*/); - geod_init(PIN->geod, PIN->a, PIN->es / (1 + sqrt (PIN->one_es))); + geod_init(PIN->geod, PIN->a, PIN->f); /* Projection specific initialization */ err = proj_errno_reset (PIN); diff --git a/src/tests/geodsigntest.c b/src/tests/geodsigntest.c index 3677e1b6e8e1..ef34b274bd0c 100644 --- a/src/tests/geodsigntest.c +++ b/src/tests/geodsigntest.c @@ -30,6 +30,19 @@ typedef double T; #define nullptr 0 #endif +#if !defined(OLD_BUGGY_REMQUO) +/* + * glibc prior to version 2.22 had a bug in remquo. This was reported in 2014 + * and fixed in 2015. See + * https://sourceware.org/bugzilla/show_bug.cgi?id=17569 + * + * The bug causes some of the tests here to fail. The failures aren't terribly + * serious (just a loss of accuracy). If you're still using the buggy glibc, + * then define OLD_BUGGY_REMQUO to be 1. + */ +#define OLD_BUGGY_REMQUO 0 +#endif + static const T wgs84_a = 6378137, wgs84_f = 1/298.257223563; /* WGS84 */ static int equiv(T x, T y) { @@ -123,7 +136,9 @@ int main() { check( geod_AngRound( 90.0 ), 90 ); checksincosd(- inf, nan, nan); +#if !OLD_BUGGY_REMQUO checksincosd(-810.0, -1.0, +0.0); +#endif checksincosd(-720.0, -0.0, +1.0); checksincosd(-630.0, +1.0, +0.0); checksincosd(-540.0, -0.0, -1.0); @@ -142,10 +157,13 @@ int main() { checksincosd(+540.0, +0.0, -1.0); checksincosd(+630.0, -1.0, +0.0); checksincosd(+720.0, +0.0, +1.0); +#if !OLD_BUGGY_REMQUO checksincosd(+810.0, +1.0, +0.0); +#endif checksincosd(+ inf, nan, nan); checksincosd( nan, nan, nan); +#if !OLD_BUGGY_REMQUO { T s1, c1, s2, c2, s3, c3; geod_sincosd( 9.0, &s1, &c1); @@ -156,6 +174,7 @@ int main() { ++n; } } +#endif check( geod_atan2d(+0.0 , -0.0 ), +180 ); check( geod_atan2d(-0.0 , -0.0 ), -180 );