From 5e7d58ec69764f7ffa8d09991aa0b0102d2f02ab Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 27 Sep 2020 15:58:54 +0200 Subject: Ortho ellipsoidal inverse: improve accuracy in polar case with (x,y) close to (0,0) --- src/projections/ortho.cpp | 4 ++-- test/gie/builtins.gie | 12 ++++++++---- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp index b4ecff7f..555f3f71 100644 --- 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 index 71622c7a..02e6766e 100644 --- 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 -- cgit v1.2.3