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.
OSGeo · Jan 24, 2025 · 560819d · 560819d
1 parent 8268efc
commit 560819d
Show file tree

Hide file tree

Showing 6 changed files with 21 additions and 72 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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