Merge pull request #2397 from cffk/merc-update

Update Mercator projection, more accurate, faster

4 files changed, 29 insertions, 48 deletions

diff --git a/src/projections/gstmerc.cpp b/src/projections/gstmerc.cpp
index 808d9ef7..50814bb5 100644
--- a/src/projections/gstmerc.cpp
+++ b/src/projections/gstmerc.cpp
@@ -28,9 +28,9 @@ static PJ_XY gstmerc_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forw
     double L, Ls, sinLs1, Ls1;
 
     L = Q->n1*lp.lam;
-    Ls = Q->c + Q->n1 * log(pj_tsfn(-1.0 * lp.phi, -1.0 * sin(lp.phi), P->e));
+    Ls = Q->c + Q->n1 * log(pj_tsfn(-lp.phi, -sin(lp.phi), P->e));
     sinLs1 = sin(L) / cosh(Ls);
-    Ls1 = log(pj_tsfn(-1.0 * asin(sinLs1), 0.0, 0.0));
+    Ls1 = log(pj_tsfn(-asin(sinLs1), -sinLs1, 0.0));
     xy.x = (Q->XS + Q->n2*Ls1) * P->ra;
     xy.y = (Q->YS + Q->n2*atan(sinh(Ls) / cos(L))) * P->ra;
 
@@ -45,9 +45,9 @@ static PJ_LP gstmerc_s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inve
 
     L = atan(sinh((xy.x * P->a - Q->XS) / Q->n2) / cos((xy.y * P->a - Q->YS) / Q->n2));
     sinC = sin((xy.y * P->a - Q->YS) / Q->n2) / cosh((xy.x * P->a - Q->XS) / Q->n2);
-    LC = log(pj_tsfn(-1.0 * asin(sinC), 0.0, 0.0));
+    LC = log(pj_tsfn(-asin(sinC), -sinC, 0.0));
     lp.lam = L / Q->n1;
-    lp.phi = -1.0 * pj_phi2(P->ctx, exp((LC - Q->c) / Q->n1), P->e);
+    lp.phi = -pj_phi2(P->ctx, exp((LC - Q->c) / Q->n1), P->e);
 
     return lp;
 }
@@ -60,13 +60,13 @@ PJ *PROJECTION(gstmerc) {
     P->opaque = Q;
 
     Q->lamc = P->lam0;
-    Q->n1 = sqrt(1.0 + P->es * pow(cos(P->phi0), 4.0) / (1.0 - P->es));
+    Q->n1 = sqrt(1 + P->es * pow(cos(P->phi0), 4.0) / (1 - P->es));
     Q->phic = asin(sin(P->phi0) / Q->n1);
-    Q->c = log(pj_tsfn(-1.0 * Q->phic, 0.0, 0.0))
-         - Q->n1 * log(pj_tsfn(-1.0 * P->phi0, -1.0 * sin(P->phi0), P->e));
-    Q->n2 = P->k0 * P->a * sqrt(1.0 - P->es) / (1.0 - P->es * sin(P->phi0) * sin(P->phi0));
+    Q->c = log(pj_tsfn(-Q->phic, -sin(P->phi0) / Q->n1, 0.0))
+         - Q->n1 * log(pj_tsfn(-P->phi0, -sin(P->phi0), P->e));
+    Q->n2 = P->k0 * P->a * sqrt(1 - P->es) / (1 - P->es * sin(P->phi0) * sin(P->phi0));
     Q->XS = 0;
-    Q->YS = -1.0 * Q->n2 * Q->phic;
+    Q->YS = -Q->n2 * Q->phic;
 
     P->inv = gstmerc_s_inverse;
     P->fwd = gstmerc_s_forward;
diff --git a/src/projections/lcc.cpp b/src/projections/lcc.cpp
index 91ffc511..46378ce4 100644
--- a/src/projections/lcc.cpp
+++ b/src/projections/lcc.cpp
@@ -106,10 +106,10 @@ PJ *PROJECTION(lcc) {
         double ml1, m1;
 
         m1 = pj_msfn(sinphi, cosphi, P->es);
-        ml1 = pj_tsfn(Q->phi1, sinphi, P->e);
-        if( ml1 == 0 ) {
+        if( fabs(Q->phi1) == M_HALFPI ) {
             return pj_default_destructor(P, PJD_ERR_LAT_1_OR_2_ZERO_OR_90);
         }
+        ml1 = pj_tsfn(Q->phi1, sinphi, P->e);
         if (secant) { /* secant cone */
             sinphi = sin(Q->phi2);
             Q->n = log(m1 / pj_msfn(sinphi, cos(Q->phi2), P->es));
@@ -117,10 +117,10 @@ PJ *PROJECTION(lcc) {
                 // Not quite, but es is very close to 1...
                 return pj_default_destructor(P, PJD_ERR_INVALID_ECCENTRICITY);
             }
-            const double ml2 = pj_tsfn(Q->phi2, sinphi, P->e);
-            if( ml2 == 0 ) {
-                return pj_default_destructor(P, PJD_ERR_LAT_1_OR_2_ZERO_OR_90);
+            if( fabs(Q->phi2) == M_HALFPI ) {
+              return pj_default_destructor(P, PJD_ERR_LAT_1_OR_2_ZERO_OR_90);
             }
+            const double ml2 = pj_tsfn(Q->phi2, sinphi, P->e);
             const double denom = log(ml1 / ml2);
             if( denom == 0 ) {
                 // Not quite, but es is very close to 1...
diff --git a/src/projections/merc.cpp b/src/projections/merc.cpp
index a77d7517..3a0ed7b4 100644
--- a/src/projections/merc.cpp
+++ b/src/projections/merc.cpp
@@ -10,45 +10,29 @@
 PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts=";
 PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Ell\n\t";
 
-#define EPS10 1.e-10
-static double logtanpfpim1(double x) {       /* log(tan(x/2 + M_FORTPI)) */
-    if (fabs(x) <= DBL_EPSILON) {
-        /* tan(M_FORTPI + .5 * x) can be approximated by  1.0 + x */
-        return log1p(x);
-    }
-    return log(tan(M_FORTPI + .5 * x));
-}
-
 static PJ_XY merc_e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
     PJ_XY xy = {0.0,0.0};
-    if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
-        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
-        return xy;
-    }
     xy.x = P->k0 * lp.lam;
-    xy.y = - P->k0 * log(pj_tsfn(lp.phi, sin(lp.phi), P->e));
+    // Instead of calling tan and sin, call sin and cos which the compiler
+    // optimizes to a single call to sincos.
+    double sphi = sin(lp.phi);
+    double cphi = cos(lp.phi);
+    xy.y = P->k0 * (asinh(sphi/cphi) - P->e * atanh(P->e * sphi));
     return xy;
 }
 
 
 static PJ_XY merc_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
     PJ_XY xy = {0.0,0.0};
-    if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
-        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
-        return xy;
-}
     xy.x = P->k0 * lp.lam;
-    xy.y = P->k0 * logtanpfpim1(lp.phi);
+    xy.y = P->k0 * asinh(tan(lp.phi));
     return xy;
 }
 
 
 static PJ_LP merc_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
     PJ_LP lp = {0.0,0.0};
-    if ((lp.phi = pj_phi2(P->ctx, exp(- xy.y / P->k0), P->e)) == HUGE_VAL) {
-        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
-        return lp;
-}
+    lp.phi = atan(pj_sinhpsi2tanphi(P->ctx, sinh(xy.y / P->k0), P->e));
     lp.lam = xy.x / P->k0;
     return lp;
 }
diff --git a/src/projections/tobmerc.cpp b/src/projections/tobmerc.cpp
index a1616036..f05a9b6b 100644
--- a/src/projections/tobmerc.cpp
+++ b/src/projections/tobmerc.cpp
@@ -9,27 +9,24 @@
 
 PROJ_HEAD(tobmerc, "Tobler-Mercator") "\n\tCyl, Sph";
 
-#define EPS10 1.e-10
-static double logtanpfpim1(double x) {       /* log(tan(x/2 + M_FORTPI)) */
-    if (fabs(x) <= DBL_EPSILON) {
-        /* tan(M_FORTPI + .5 * x) can be approximated by  1.0 + x */
-        return log1p(x);
-    }
-    return log(tan(M_FORTPI + .5 * x));
-}
-
 static PJ_XY tobmerc_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
     PJ_XY xy = {0.0, 0.0};
     double cosphi;
 
-    if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
+    if (fabs(lp.phi) >= M_HALFPI) {
+      // builtins.gie tests "Test expected failure at the poles:".  However
+      // given that M_HALFPI is strictly less than pi/2 in double precision,
+      // it's not clear why shouldn't just return a large result for xy.y (and
+      // it's not even that large, merely 38.025...).  Even if the logic was
+      // such that phi was strictly equal to pi/2, allowing xy.y = inf would be
+      // a reasonable result.
         proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
         return xy;
     }
 
     cosphi = cos(lp.phi);
     xy.x = P->k0 * lp.lam * cosphi * cosphi;
-    xy.y = P->k0 * logtanpfpim1(lp.phi);
+    xy.y = P->k0 * asinh(tan(lp.phi));
     return xy;
 }