Merge pull request #2397 from cffk/merc-update

Update Mercator projection, more accurate, faster

1 files changed, 23 insertions, 9 deletions

diff --git a/src/tsfn.cpp b/src/tsfn.cpp
index 32da09f2..8ed258d6 100644
--- a/src/tsfn.cpp
+++ b/src/tsfn.cpp
@@ -4,14 +4,28 @@
 #include "proj_internal.h"
 
 double pj_tsfn(double phi, double sinphi, double e) {
-    double denominator;
-    sinphi *= e;
+  /****************************************************************************
+   * Determine function ts(phi) defined in Snyder (1987), Eq. (7-10)
+   * Inputs:
+   *   phi = geographic latitude (radians)
+   *   e = eccentricity of the ellipsoid (dimensionless)
+   * Output:
+   *   ts = exp(-psi) where psi is the isometric latitude (dimensionless)
+   *      = 1 / (tan(chi) + sec(chi))
+   * Here isometric latitude is defined by
+   *   psi = log( tan(pi/4 + phi/2) *
+   *              ( (1 - e*sin(phi)) / (1 + e*sin(phi)) )^(e/2) )
+   *       = asinh(tan(phi)) - e * atanh(e * sin(phi))
+   *       = asinh(tan(chi))
+   *   chi = conformal latitude
+   ***************************************************************************/
 
-    /* avoid zero division, fail gracefully */
-    denominator = 1.0 + sinphi;
-    if (denominator == 0.0)
-        return HUGE_VAL;
-
-    return (tan (.5 * (M_HALFPI - phi)) /
-            pow((1. - sinphi) / (denominator), .5 * e));
+    double cosphi = cos(phi);
+    // exp(-asinh(tan(phi))) = 1 / (tan(phi) + sec(phi))
+    //                       = cos(phi) / (1 + sin(phi)) good for phi > 0
+    //                       = (1 - sin(phi)) / cos(phi) good for phi < 0
+    return exp(e * atanh(e * sinphi)) *
+      ( sinphi > 0 ?
+        cosphi / (1 + sinphi) :
+        (1 - sinphi) / cosphi );
 }