Robinson: fix wrong values for forward path for latitudes >= 87.5 (fixes

#1172), and fix inaccurate inverse method
OSGeo · Feb 24, 2019 · 9118523 · 9118523
1 parent 38a1525
commit 9118523
Show file tree

Hide file tree

Showing 2 changed files with 29 additions and 7 deletions.
diff --git a/src/projections/robin.cpp b/src/projections/robin.cpp
@@ -7,7 +7,7 @@
 PROJ_HEAD(robin, "Robinson") "\n\tPCyl, Sph";
 
 #define V(C,z) (C.c0 + z * (C.c1 + z * (C.c2 + z * C.c3)))
-#define DV(C,z) (C.c1 + z * (C.c2 + C.c2 + z * 3. * C.c3))
+#define DV(C,z) (C.c1 + 2 * z * C.c2 + z * z * 3. * C.c3)
 
 /*
 note: following terms based upon 5 deg. intervals in degrees.
@@ -74,7 +74,7 @@ static const struct COEFS Y[] = {
 #define RC1 0.08726646259971647884
 #define NODES   18
 #define ONEEPS  1.000001
-#define EPS 1e-8
+#define EPS 1e-10
 /* Not sure at all of the appropriate number for MAX_ITER... */
 #define MAX_ITER 100
 
@@ -90,7 +90,7 @@ static PJ_XY s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
         proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
         return xy;
     }
-    if (i >= NODES) i = NODES - 1;
+    if (i >= NODES) i = NODES;
     dphi = RAD_TO_DEG * (dphi - RC1 * i);
     xy.x = V(X[i], dphi) * FXC * lp.lam;
     xy.y = V(Y[i], dphi) * FYC;
@@ -133,10 +133,8 @@ static PJ_LP s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inverse */
         T = Y[i];
         /* first guess, linear interp */
         t = 5. * (lp.phi - T.c0)/(Y[i+1].c0 - T.c0);
-        /* make into root */
-        T.c0 = (float)(T.c0 - lp.phi);
         for (iters = MAX_ITER; iters ; --iters) { /* Newton-Raphson */
-            t -= t1 = V(T,t) / DV(T,t);
+            t -= t1 = (V(T,t) - lp.phi) / DV(T,t);
             if (fabs(t1) < EPS)
                 break;
         }

diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
@@ -4438,7 +4438,7 @@ Robinson
 ===============================================================================
 
 -------------------------------------------------------------------------------
-operation +proj=robin   +a=6400000    +lat_1=0.5 +lat_2=2
+operation +proj=robin   +a=6400000
 -------------------------------------------------------------------------------
 tolerance 0.1 mm
 accept  2 1
@@ -4450,6 +4450,18 @@ expect  -189588.423282508 107318.530350703
 accept  -2 -1
 expect  -189588.423282508 -107318.530350703
 
+accept  0 89.5
+expect  0.000000000000   8639799.718722090125
+
+accept  0 90
+expect  0.000000000000   8654720.000000000000
+
+accept  0 -89.5
+expect  0.000000000000   -8639799.718722090125
+
+accept  0 -90
+expect  0.000000000000   -8654720.000000000000
+
 direction inverse
 accept  200 100
 expect  0.002109689 0.000931806
@@ -4460,6 +4472,18 @@ expect  -0.002109689 0.000931806
 accept  -200 -100
 expect  -0.002109689 -0.000931806
 
+accept  0.000000000000   8639799.718722090125
+expect  0 89.5
+
+accept  0.000000000000   8654720.000000000000
+expect  0 90
+
+accept  0.000000000000   -8639799.718722090125
+expect  0 -89.5
+
+accept  0.000000000000   -8654720.000000000000
+expect  0 -90
+
 
 ===============================================================================
 Roussilhe Stereographic