Ortho ellipsoidal inverse: improve accuracy in polar case with (x,y) …

…close to (0,0)
OSGeo · Sep 27, 2020 · 5e7d58e · 5e7d58e
1 parent fedeeec
commit 5e7d58e
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 6 deletions.
diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp
@@ -171,7 +171,7 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) {           /* Ellipsoidal, inver
         // ==> lam = atan2(xy.x, -xy.y * sign(phi0))
         // ==> xy.x^2 + xy.y^2 = nu^2 * cosphi^2
         //                rh^2 = cosphi^2 / (1 - es * sinphi^2)
-        // ==>  sinphi^2 = (1 - rh^2) / (1 - es * rh^2)
+        // ==>  cosphi^2 = rh^2 * (1 - es) / (1 - es * rh^2)
 
         const double rh2 = SQ(xy.x) + SQ(xy.y);
         if (rh2 >= 1. - 1e-15) {
@@ -184,7 +184,7 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) {           /* Ellipsoidal, inver
             lp.phi = 0;
         }
         else {
-            lp.phi = asin(sqrt((1 - rh2) / (1 - P->es * rh2))) * (Q->mode == N_POLE ? 1 : -1);
+            lp.phi = acos(sqrt(rh2 * P->one_es / (1 - P->es * rh2))) * (Q->mode == N_POLE ? 1 : -1);
         }
         lp.lam = atan2(xy.x, xy.y * (Q->mode == N_POLE ? -1 : 1));
         return lp;

diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
@@ -4565,8 +4565,12 @@ accept      0               90
 expect      0               0
 roundtrip   1
 
-accept      1               89.999
-expect      1.9493         -111.6770
+accept      30              45
+expect      2258795.4394   -3912348.4650
+roundtrip   1
+
+accept      135             89.999999873385
+expect      0.01            0.01
 roundtrip   1
 
 # Point not visible from the projection plane
@@ -4592,8 +4596,8 @@ accept      0               -90
 expect      0               0
 roundtrip   1
 
-accept      1               -89.999
-expect      1.9493          111.6770
+accept      135             -89.999999873385
+expect      0.01            -0.01
 roundtrip   1
 
 # Point not visible from the projection plane