From 560819d4cd4fb32d515d961df86e8b2e8ceb2c2c Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Fri, 24 Jan 2025 14:21:37 -0500 Subject: [PATCH] Use C++17 function ellint_1 in adams.cpp Replace approximate series for F(x, 1/sqrt(2)) by a call to ellint_1. The approximate series was only accurate to "better than 1e-7", which meant that the errors in the projection were typically about 1 m. The adoption of ellint_1 (with its considerably better accuracy) then leads to a torrent of errors in the tests. These have been fixed by changing the tolerances from 1 mm to 1 m. The values in the tests should be replaced by more accurate ones, but I'm reluctant merely to update the values with the ones the code now produces given that the code depends on a parameter TOL = 1e-9. Also, I don't understand why K(1/sqrt(2)) = std::comp_ellint_1(RSQRT2) = 1.8540746773013719 is sometimes given to only 5 decimal places as 1.85407. --- src/projections/adams.cpp | 37 ++++++-------------------------- test/gie/adams_hemi.gie | 2 +- test/gie/adams_ws1.gie | 2 +- test/gie/adams_ws2.gie | 6 ++---- test/gie/guyou.gie | 2 +- test/gie/peirce_q.gie | 44 +++++++++------------------------------ 6 files changed, 21 insertions(+), 72 deletions(-) diff --git a/src/projections/adams.cpp b/src/projections/adams.cpp index cdef83ec46..9f0845a146 100644 --- a/src/projections/adams.cpp +++ b/src/projections/adams.cpp @@ -78,34 +78,7 @@ struct pj_adams_data { } // anonymous namespace #define TOL 1e-9 -#define RSQRT2 0.7071067811865475244008443620 - -static double ell_int_5(double phi) { - /* Procedure to compute elliptic integral of the first kind - * where k^2=0.5. Precision good to better than 1e-7 - * The approximation is performed with an even Chebyshev - * series, thus the coefficients below are the even values - * and where series evaluation must be multiplied by the argument. */ - constexpr double C0 = 2.19174570831038; - static const double C[] = { - -8.58691003636495e-07, 2.02692115653689e-07, 3.12960480765314e-05, - 5.30394739921063e-05, -0.0012804644680613, -0.00575574836830288, - 0.0914203033408211, - }; - - double y = phi * M_2_PI; - y = 2. * y * y - 1.; - double y2 = 2. * y; - double d1 = 0.0; - double d2 = 0.0; - for (double c : C) { - double temp = d1; - d1 = y2 * d1 - d2 + c; - d2 = temp; - } - - return phi * (y * d1 - d2 + 0.5 * C0); -} +#define RSQRT2 0.7071067811865475244008443621 static PJ_XY adams_forward(PJ_LP lp, PJ *P) { double a = 0., b = 0.; @@ -198,14 +171,16 @@ static PJ_XY adams_forward(PJ_LP lp, PJ *P) { if (sn) n = -n; - xy.x = ell_int_5(m); - xy.y = ell_int_5(n); + // Elliptic integral of the first kind with k^2 = 0.5 + xy.x = std::ellint_1(RSQRT2, m); + xy.y = std::ellint_1(RSQRT2, n); if (Q->mode == PEIRCE_Q) { /* Constant complete elliptic integral of the first kind with m=0.5, * calculated using * https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ellipk.html - * . Used as basic as scaled shift distance */ + * . Used as basic as scaled shift distance. + * Could replace 1.85407... by std::comp_ellint_1(RSQRT2) */ constexpr double shd = 1.8540746773013719 * 2; /* For square and diamond Quincuncial projections, spin out southern diff --git a/test/gie/adams_hemi.gie b/test/gie/adams_hemi.gie index 9d23b939a3..14cfd163a3 100644 --- a/test/gie/adams_hemi.gie +++ b/test/gie/adams_hemi.gie @@ -6,7 +6,7 @@ ------------------------------------------------------------ operation +proj=adams_hemi +R=6370997 -tolerance 1 mm +tolerance 1 m ------------------------------------------------------------ accept -179.0512914938 -90.1445918836 expect failure errno coord_transfm_invalid_coord diff --git a/test/gie/adams_ws1.gie b/test/gie/adams_ws1.gie index 3593fe12e1..636ee9a3f3 100644 --- a/test/gie/adams_ws1.gie +++ b/test/gie/adams_ws1.gie @@ -7,7 +7,7 @@ ------------------------------------------------------------ operation +proj=adams_ws1 +R=6370997 -tolerance 1 mm +tolerance 1 m ------------------------------------------------------------ accept -179.5170670673 -90.3642618405 expect failure errno coord_transfm_invalid_coord diff --git a/test/gie/adams_ws2.gie b/test/gie/adams_ws2.gie index 8ba0475559..6c02b5534b 100644 --- a/test/gie/adams_ws2.gie +++ b/test/gie/adams_ws2.gie @@ -8,7 +8,7 @@ ------------------------------------------------------------ operation +proj=adams_ws2 +R=6370997 -tolerance 1 mm +tolerance 1 m ------------------------------------------------------------ accept -179.7092450238 -90.0290393775 expect failure errno coord_transfm_invalid_coord @@ -2123,9 +2123,9 @@ expect failure errno coord_transfm_invalid_coord # Test inverse ------------------------------------------------------------------------------- operation +proj=adams_ws2 +ellps=WGS84 +tolerance 1 m ------------------------------------------------------------------------------- direction forward -tolerance 1 mm accept 0 0 expect 0 0 @@ -2151,8 +2151,6 @@ accept 0 -89.999 expect 0 -15743336.122 roundtrip 1 -# Results a bit different on x86 -tolerance 3 mm accept 179.999 89.999 expect 693320.704 16030515.906 roundtrip 1 diff --git a/test/gie/guyou.gie b/test/gie/guyou.gie index 2ed2ab49b8..56e716d976 100644 --- a/test/gie/guyou.gie +++ b/test/gie/guyou.gie @@ -6,7 +6,7 @@ ------------------------------------------------------------ operation +proj=guyou +R=6370997 -tolerance 1 mm +tolerance 1 m ------------------------------------------------------------ accept -179.2338274749 -90.7265739758 expect failure errno coord_transfm_invalid_coord diff --git a/test/gie/peirce_q.gie b/test/gie/peirce_q.gie index 0fac30edcf..633680a1a2 100644 --- a/test/gie/peirce_q.gie +++ b/test/gie/peirce_q.gie @@ -15,7 +15,7 @@ ------------------------------------------------------------ operation +proj=peirce_q +R=6370997 +shape=square -tolerance 10 mm +tolerance 1 m ------------------------------------------------------------ accept -179.6126302052 -90.2440064745 @@ -139,7 +139,7 @@ expect -153172.25 215556.52 ------------------------------------------------------------ operation +proj=peirce_q +R=6370997 +shape=diamond - tolerance 10 mm + tolerance 1 m ------------------------------------------------------------ accept -179.6126302052 -90.2440064745 @@ -263,7 +263,7 @@ expect 44112.34 260730.62 ------------------------------------------------------------ operation +proj=peirce_q +R=6370997 +shape=horizontal -tolerance 10 mm +tolerance 1 m ------------------------------------------------------------ accept -179.6126302052 -90.2440064745 @@ -387,7 +387,7 @@ expect -11768191.87 260730.62 ------------------------------------------------------------ operation +proj=peirce_q +R=6370997 +shape=horizontal +scrollx=0.75 - tolerance 10 mm + tolerance 1 m ------------------------------------------------------------ accept -179.6126302052 -90.2440064745 @@ -511,7 +511,7 @@ expect -23580496.07 260730.62 ------------------------------------------------------------ operation +proj=peirce_q +R=6370997 +shape=vertical - tolerance 10 mm + tolerance 1 m ------------------------------------------------------------ accept -179.6126302052 -90.2440064745 @@ -635,7 +635,7 @@ expect 44112.34 -11551573.59 ------------------------------------------------------------ operation +proj=peirce_q +R=6370997 +shape=vertical +scrolly=-0.25 - tolerance 10 mm + tolerance 1 m ------------------------------------------------------------ accept -179.6126302052 -90.2440064745 @@ -759,7 +759,7 @@ expect 44112.34 -23363877.80 ------------------------------------------------------------ operation +proj=peirce_q +R=6370997 +shape=nhemisphere - tolerance 10 mm + tolerance 1 m ------------------------------------------------------------ accept -179.6126302052 -90.2440064745 @@ -883,7 +883,7 @@ expect 44112.34 260730.62 ------------------------------------------------------------ operation +proj=peirce_q +R=6370997 +shape=shemisphere - tolerance 10 mm + tolerance 1 m ------------------------------------------------------------ accept -179.6126302052 -90.2440064745 @@ -1008,12 +1008,9 @@ expect failure errno coord_transfm_outside_projection_domain # Test inverse ------------------------------------------------------------ operation +proj=peirce_q +shape=square +tolerance 1 m ------------------------------------------------------------ -#tolerance 1 mm -# has to bump to this for i386 -tolerance 150 mm - accept 0 90 expect 0 0 roundtrip 1 @@ -1024,9 +1021,7 @@ roundtrip 1 accept 0 -90 expect 16723842.303160080686 -16723842.303160080686 -#tolerance 2 mm roundtrip 1 -#tolerance 1 mm accept 0 45 expect 3725360.212758612353 -3725360.212758612353 @@ -1037,10 +1032,8 @@ expect 12998482.090401465073 -12998482.090401465073 roundtrip 1 accept 45 0 -tolerance 200 mm expect 16723842.564696932212 -0.095041956369 roundtrip 1 -tolerance 150 mm accept -45 0 expect 0 -16723842.469654975459 @@ -1169,21 +1162,18 @@ expect -16723763.588384689763 16723763.588384689763 # Test inverse ------------------------------------------------------------ operation +proj=peirce_q +shape=diamond +tolerance 1 m ------------------------------------------------------------ -#tolerance 1 mm # has to bump to this for i386 -tolerance 150 mm accept 0 90 expect 0 0 roundtrip 1 accept 0 -90 -#tolerance 10 mm expect 0 -23651084.600117880851 roundtrip 1 -#tolerance 1 mm accept 0 45 expect 0 -5268454.937608348206 @@ -1194,10 +1184,8 @@ expect 0 -18382629.662509534508 roundtrip 1 accept 45 0 -tolerance 200 mm expect 11825542.417788611725 -11825542.552198234946 roundtrip 1 -tolerance 150 mm accept -45 0 expect -11825542.417788611725 -11825542.417788611725 @@ -1209,9 +1197,6 @@ roundtrip 1 accept -90 0 expect -11825542.417788611725 0.000000000000 -#tolerance 20 mm -#roundtrip 1 -#tolerance 1 mm accept 135 0 expect 11825542.552198234946 11825542.417788611725 @@ -1223,15 +1208,9 @@ roundtrip 1 accept 179.99 0 expect 1574.295465656175 11825542.417788611725 -#tolerance 200 mm -#roundtrip 1 -#tolerance 1 mm accept -179.99 0 expect -1574.295465656175 11825542.417788611725 -#tolerance 30 mm -#roundtrip 1 -#tolerance 1 mm accept 45 45 expect 3747362.066324858926 -3747362.066324859392 @@ -1279,9 +1258,6 @@ roundtrip 1 accept -90 -45 expect -18382629.662509530783 0.000000000000 -#tolerance 3 mm -#roundtrip 1 -#tolerance 1 mm accept 135 -45 expect 3747362.066324858926 19903722.533793020993