Ortho ellipsoidal inverse: add non iterative implementations for polar and equatorial

2 files changed, 111 insertions, 8 deletions

diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp
index 0594205b..ac54d88a 100644
--- a/src/projections/ortho.cpp
+++ b/src/projections/ortho.cpp
@@ -160,6 +160,66 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) {           /* Ellipsoidal, inver
     PJ_LP lp;
     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
 
+    if( Q->mode == N_POLE || Q->mode == S_POLE )
+    {
+        // Polar case. Forward case equations can be simplified as:
+        // xy.x = nu * cosphi * sinlam
+        // xy.y = nu * -cosphi * coslam * sign(phi0)
+        // ==> 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)
+
+        const double rh2 = xy.x * xy.x + xy.y * xy.y;
+        if (rh2 >= 1. - 1e-15) {
+            if ((rh2 - 1.) > EPS10) {
+                proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+                proj_log_trace(P, "Point (%.3f, %.3f) is outside the projection boundary");
+                lp.lam = HUGE_VAL; lp.phi = HUGE_VAL;
+                return lp;
+            }
+            lp.phi = 0;
+        }
+        else {
+            lp.phi = asin(sqrt((1 - rh2) / (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;
+    }
+
+    if( Q->mode == EQUIT )
+    {
+        // Equatorial case. Forward case equations can be simplified as:
+        // xy.x = nu * cosphi * sinlam
+        // xy.y  = nu * sinphi * (1 - P->es)
+        // x^2 * (1 - es * sinphi^2) = (1 - sinphi^2) * sinlam^2
+        // y^2 / ((1 - es)^2 + y^2 * es) = sinphi^2
+
+        const auto SQ = [](double x) { return x*x; };
+
+        // Equation of the ellipse
+        if( SQ(xy.x) + SQ(xy.y * (P->a / P->b)) > 1 + 1e-11 ) {
+            proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+            proj_log_trace(P, "Point (%.3f, %.3f) is outside the projection boundary");
+            lp.lam = HUGE_VAL; lp.phi = HUGE_VAL;
+            return lp;
+        }
+
+        const double sinphi2 = xy.y == 0 ? 0 : 1.0 / (SQ((1 - P->es) / xy.y) + P->es);
+        if( sinphi2 > 1 - 1e-11 ) {
+            lp.phi = M_PI_2 * (xy.y > 0 ? 1 : -1);
+            lp.lam = 0;
+            return lp;
+        }
+        lp.phi = asin(sqrt(sinphi2)) * (xy.y > 0 ? 1 : -1);
+        const double sinlam = xy.x * sqrt((1 - P->es * sinphi2) / (1 - sinphi2));
+        if( fabs(sinlam) - 1 > -1e-15 )
+            lp.lam = M_PI_2 * (xy.x > 0 ? 1: -1);
+        else
+            lp.lam = asin(sinlam);
+        return lp;
+    }
+
     // From EPSG guidance note 7.2
 
     // It suggests as initial guess:
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index 9d5b6644..b63b7902 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -4270,7 +4270,7 @@ expect  -223374.577355253 -111701.072127637
 ===============================================================================
 
 -------------------------------------------------------------------------------
-# Test the equatorial aspect of the Orthopgraphic projection.
+# Test the equatorial aspect of the Orthographic projection.
 
 # Test data from Snyder (1987), table 22, p. 151.
 -------------------------------------------------------------------------------
@@ -4322,7 +4322,7 @@ expect      failure errno tolerance_condition
 
 
 -------------------------------------------------------------------------------
-# Test the oblique aspect of the Orthopgraphic projection.
+# Test the oblique aspect of the Orthographic projection.
 
 # Test data from Snyder (1987), table 23, pp. 152-153.
 -------------------------------------------------------------------------------
@@ -4356,7 +4356,7 @@ expect      failure errno tolerance_condition
 
 
 -------------------------------------------------------------------------------
-# Test the north polar aspect of the Orthopgraphic projection.
+# Test the north polar aspect of the Orthographic projection.
 -------------------------------------------------------------------------------
 operation   +proj=ortho +R=1 +lat_0=90 +lon_0=0
 -------------------------------------------------------------------------------
@@ -4386,7 +4386,7 @@ accept      2   2
 expect      failure errno tolerance_condition
 
 -------------------------------------------------------------------------------
-# Test the south polar aspect of the Orthopgraphic projection.
+# Test the south polar aspect of the Orthographic projection.
 -------------------------------------------------------------------------------
 operation   +proj=ortho +R=1 +lat_0=-90 +lon_0=0
 -------------------------------------------------------------------------------
@@ -4443,6 +4443,22 @@ accept      0               0
 expect      0               0
 roundtrip   1
 
+accept      1               1
+expect      111296.9991     110568.7748
+roundtrip   1
+
+accept      1               -1
+expect      111296.9991     -110568.7748
+roundtrip   1
+
+accept      -1              1
+expect      -111296.9991    110568.7748
+roundtrip   1
+
+accept      -1              -1
+expect      -111296.9991    -110568.7748
+roundtrip   1
+
 accept      89.99           0
 expect      6378136.9029    0
 roundtrip   1
@@ -4468,20 +4484,37 @@ accept      -90.00001       0
 expect      failure         errno tolerance_condition
 
 # Consistant with WGS84 semi-major axis
-# The inverse transformation doesn't converge due to properties of the projection
 accept      90              0
 expect      6378137         0
+roundtrip   1
 
 accept      -90             0
 expect      -6378137        0
+roundtrip   1
 
 # Consistant with WGS84 semi-minor axis
-# The inverse transformation doesn't converge due to properties of the projection
 accept      0               90
 expect      0               6356752.3142
+roundtrip   1
 
 accept      0               -90
 expect      0               -6356752.3142
+roundtrip   1
+
+# Point not visible from the projection plane
+direction   inverse
+accept      0               6356752.3143
+expect      failure         errno tolerance_condition
+
+# Point not visible from the projection plane
+direction   inverse
+accept      1000            6356752.314
+expect      failure         errno tolerance_condition
+
+# Point not visible from the projection plane
+direction   inverse
+accept      6378137.0001    0
+expect      failure         errno tolerance_condition
 
 -------------------------------------------------------------------------------
 # North pole tests
@@ -4501,9 +4534,14 @@ accept      0               -0.0000001
 expect      failure         errno tolerance_condition
 
 # Consistant with WGS84 semi-major axis
-# The inverse transformation doesn't converge due to properties of the projection
 accept      0               0
 expect      0               -6378137
+roundtrip   1
+
+# Point not visible from the projection plane
+direction   inverse
+accept      0               -6378137.1
+expect      failure         errno tolerance_condition
 
 -------------------------------------------------------------------------------
 # South pole tests
@@ -4523,9 +4561,14 @@ accept      0               0.0000001
 expect      failure         errno tolerance_condition
 
 # Consistant with WGS84 semi-major axis
-# The inverse transformation doesn't converge due to properties of the projection
 accept      0               0
 expect      0               6378137
+roundtrip   1
+
+# Point not visible from the projection plane
+direction   inverse
+accept      0               6378137.1
+expect      failure         errno tolerance_condition
 
 
 ===============================================================================