cass: fix forward computation of easting (fixes OSGeo#3432), and use …

…iterative generic inversion to make GIGS 5108 roundtripping tests succeed

The error on the sign after "nu * A * (1. - A2 * T * (C1 ...") dates
back to PROJ origins. All other independent implementations (Apache SIS
cf
/~https://github.com/apache/sis/blob/a04bd4b29816c91fc4dfc66135b8de5e962933a0/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CassiniSoldner.java,
CSMAP cf
/~https://github.com/normanb/CS-Map/blob/1d8f2db7f13194a99dc218d4cd91754c64a2a774/Source/CS_csini.c)
use the same formulas as EPSG guidance note which is consistent with
Snyder's Map projections - A working manual). So align on them, so that
everyone is consistent (and hopefully right).

While doing it, I've modified the inverse method to use the generic
inversion method to improve its accuracy, which make it possible to
enable the 5108.gie tests that were previously failing on the 1000
iterations roundtripping.

src/generic_inverse.cpp

@@ -51,7 +51,7 @@ PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial) {
         PJ_XY xyApprox = P->fwd(lp, P);
         const double deltaX = xyApprox.x - xy.x;
         const double deltaY = xyApprox.y - xy.y;
-        if (fabs(deltaX) < 1e-10 && fabs(deltaY) < 1e-10) {
+        if (fabs(deltaX) < 1e-12 && fabs(deltaY) < 1e-12) {
             return lp;
         }
 

src/projections/cass.cpp

@@ -43,7 +43,7 @@ static PJ_XY cass_e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward
     const double A2 = A * A;
 
     xy.x = nu * A * (1. - A2 * T *
-        (C1 - (8. - T + 8. * C) * A2 * C2));
+        (C1 + (8. - T + 8. * C) * A2 * C2));
     xy.y = M - Q->m0 + nu * tanphi * A2 *
         (.5 + (5. - T + 6. * C) * A2 * C3);
     if( Q->hyperbolic )
@@ -82,13 +82,12 @@ static PJ_LP cass_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse
     lp.lam = D * (1. + T1 * D2 *
         (-C4 + (1. + 3. * T1) * D2 * C5)) / cos (phi1);
 
-    if( Q->hyperbolic )
-    {
-        // EPSG guidance note 7-2 suggests a custom approximation for the
-        // 'Vanua Levu 1915 / Vanua Levu Grid' case, but better use the
-        // generic inversion method
-        lp = pj_generic_inverse_2d(xy, P, lp);
-    }
+    // EPSG guidance note 7-2 suggests a custom approximation for the
+    // 'Vanua Levu 1915 / Vanua Levu Grid' case, but better use the
+    // generic inversion method
+    // Actually use it in the non-hyperbolic case. It enables to make the
+    // 5108.gie roundtripping tests to success, with at most 2 iterations.
+    lp = pj_generic_inverse_2d(xy, P, lp);
 
     return lp;
 }

test/CMakeLists.txt

@@ -80,7 +80,7 @@ proj_add_gie_test("GIGS-5104" "gigs/5104.gie")
 proj_add_gie_test("GIGS-5105.2" "gigs/5105.2.gie")
 proj_add_gie_test("GIGS-5106" "gigs/5106.gie")
 proj_add_gie_test("GIGS-5107" "gigs/5107.gie")
-#proj_add_gie_test("GIGS-5108" "gigs/5108.gie")
+proj_add_gie_test("GIGS-5108" "gigs/5108.gie")
 proj_add_gie_test("GIGS-5109" "gigs/5109.gie")
 #proj_add_gie_test("GIGS-5110" "gigs/5110.gie")
 proj_add_gie_test("GIGS-5111.1" "gigs/5111.1.gie")

test/gie/builtins.gie

@@ -850,6 +850,18 @@ expect  -0.001790493 0.000895247
 accept  -200 -100
 expect  -0.001790493 -0.000895247
 
+-------------------------------------------------------------------------------
+# test point from EPSG Guidance Note 7.2
+operation +proj=cass +lat_0=10.4416666666667 +lon_0=-61.3333333333333 \
+          +x_0=86501.46392052 +y_0=65379.0134283 \
+          +a=6378293.64520876 +b=6356617.98767984 +to_meter=0.201166195164
+-------------------------------------------------------------------------------
+
+tolerance 0.1 mm
+accept  -62                10
+expect  66644.94040882     82536.21873655
+roundtrip 1
+
 -------------------------------------------------------------------------------
 # Hyperbolic variant: test point from EPSG Guidance Note 7.2
 operation +proj=cass +hyperbolic +a=6378306.376305601 +rf=293.466307 \
@@ -859,7 +871,7 @@ operation +proj=cass +hyperbolic +a=6378306.376305601 +rf=293.466307 \
 
 tolerance 0.1 mm
 accept  179.99433652777776 -16.841456527777776
-expect  16015.289017555102 13369.660053668682
+expect  16015.28901692     13369.66005367
 roundtrip 1
 
 ===============================================================================

test/gigs/5108.gie.failing → test/gigs/5108.gie

@@ -4,18 +4,19 @@ Test 5108, Cassini-Soldner, v2-0_2011-06-28.
 
 --------------------------------------------------------------------------------
 
-# GDM2000
-<4742> +proj=longlat +ellps=GRS80  <>
+<gie-strict>
 
-# GDM2000 / Johor Grid
-<3377> +proj=cass +lat_0=2.121679744444445 +lon_0=103.4279362361111 +x_0=-14810.562 +y_0=8758.32 +ellps=GRS80 +units=m  <>
+use_proj4_init_rules true
 
+# GDM2000
+# <4742> +proj=longlat +ellps=GRS80  <>
 
-<gie>
+# GDM2000 / Johor Grid
+# <3377> +proj=cass +lat_0=2.121679744444445 +lon_0=103.4279362361111 +x_0=-14810.562 +y_0=8758.32 +ellps=GRS80 +units=m  <>
 
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
-           +step +init=epsg:4742 +inv
+operation  +proj=pipeline \
+           +step +init=epsg:4742 +inv \
            +step +init=epsg:3377
 --------------------------------------------------------------------------------
 tolerance  0.05 m
@@ -87,8 +88,8 @@ accept    104 5
 expect    48630.563 327067.097
 
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
-           +step +init=epsg:3377 +inv
+operation  +proj=pipeline \
+           +step +init=epsg:3377 +inv \
            +step +init=epsg:4742
 --------------------------------------------------------------------------------
 tolerance  0.05 m
@@ -160,8 +161,8 @@ accept    48630.563 327067.097
 expect    104 5
 
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
-           +step +init=epsg:4742 +inv
+operation  +proj=pipeline \
+           +step +init=epsg:4742 +inv \
            +step +init=epsg:3377
 --------------------------------------------------------------------------------
 tolerance  0.006 m
@@ -232,4 +233,4 @@ tolerance  0.006 m
 accept    104 5
 roundtrip 1000
 
-</gie>
+</gie-strict>