Update Mercator projection

Introduction
------------

The existing formulation for the Mercator projection is
"satisfactory"; it is reasonably accurate.  However for a core
projection like Mercator, I think we should strive for full double
precision accuracy.

This commit uses cleaner, more accurate, and faster methods for
computing the forward and inverse projections.  These use the
formulation in terms of hyperbolic functions that are manifestly odd
in latitude

    psi = asinh(tan(phi)) - e * atanh(e  * sin(phi))

(phi = latitude; psi = isometric latitude = Mercator y coordinate).
Contrast this with the existing formulation

    psi = log(tan(pi/4 - phi/2))
          - e/2 * log((1 + e * sin(phi)) / (1 - e * sin(phi)))

where psi(-phi) isn't exactly equal to -psi(phi) and psi(0) isn't
guaranteed to be 0.

Implementation
--------------

There's no particular issue implementing the forward projection, just
apply the formulas above.  The inverse projection is tricky because
there's no closed form solution for the inverse.  The existing code
for the inverse uses an iterative method from Snyder.  This is the
usual hokey function iteration, and, as usual, the convergence rate is
linear (error reduced by a constant factor on each iteration).  This
is OK (just) for low accuracy work.  But nowadays, something with
quadratic convergence (e.g., Newton's method, number of correct digits
doubles on each iteration) is preferred (and used here).  More on this
later.

The solution for phi(psi) I use is described in my TM paper and I
lifted the specific formulation from GeographicLib's Math::tauf, which
uses the same underlying machinery for all conformal projections.  It
solves for tan(phi) in terms of sinh(psi) which as a near identity
mapping is ideal for Newton's method.

For comparison I also look at the approach adopted by Poder + Engsager
in their TM paper and implemented in etmerc.  This uses trigonometric
series (accurate to n^6) to convert phi <-> chi.  psi is then given by

    psi = asinh(tan(chi))

Accuracy
--------

I tested just the routines for transforming phi <-> psi from merc.cpp
and measured the errors (converted to true nm = nanometers) for the
forward and inverse mapping.  I also included in my analysis the
method used by etmerc.  This uses a trigonometric series to convert
phi <-> chi = atan(sinh(psi)), the conformal latitude.

		forward      inverse
		max  rms     max    rms
    old merc    3.60 0.85 2189.47 264.81
      etmerc    1.82 0.38    1.42   0.37
    new merc    1.83 0.30    2.12   0.31

1 nm is pretty much the absolute limit for accuracy in double
precision (1 nm = 10e6 m / 2^53, approximately), and 5 nm is probably
the limit on what you should routinely expect.  So the old merc
inverse is considerably less accurate that it could be.  The old merc
forward is OK on accuracy -- except that if does not preserve the
parity of the projection.

The accuracy of etmerc is fine (the truncation error of the 6th order
series is small compared with the round-off error).  However,
situation reverses as the flattening is increased.  E.g., at f =
1/150, the max error for the inverse projection is 8 nm.  etmerc is OK
for terrestrial applications, but couldn't be used for Mars.

Timing
------

Here's what I get with g++ -O3 on various Linux machines with recent
versions of g++.  As always, you should take these with a grain of
salt.  You might expect the relative timings to vary by 20% or so when
switching between compilers/machines.  Times per call in ns =
nanoseconds.

               forward inverse
    old merc     121      360
      etmerc    4e-6      1.4
    new merc      20      346

The new merc method is 6 times faster at the forward projection and
modestly faster at the inverse projection (despite being more
accurate).  The latter result is because it only take 2 iterations of
Newton's method to get full accuracy compared with an average of 5
iterations for the old method to get only um accuracy.

A shocking aspect of these timings is how fast etmerc is.  Another is
that forward etmerc is streaks faster that inverse etmerc (it made be
doubt my timing code).  Evidently, asinh(tan(chi)) is a lot faster to
compute than atan(sinh(psi)).  The hesitation about adopting etmerc
then comes down to:

  * the likelihood that Mercator may be used for non-terrestrial
    bodies;

  * the question of whether the timing benefits for the etmerc method
    would be noticeable in a realistic application;

  * need to duplicate the machinery for evaluating the coefficients
    for the series and for Clenshaw summation in the current code
    layout.

Ripple effects
==============

The Mercator routines used the the Snyder method, pj_tsfn and pj_phi2,
are used in other projections.  These relate phi to t = exp(-psi) (a
rather bizarre choice in my book).  I've retrofitted these to use the
more accurate methods.  These do the "right thing" for phi in [-pi/2,
pi/2] , t in [0, inf], and e in [0, 1).  NANs are properly handled.

Of course, phi = pi/2 in double precision is actually less than pi/2,
so cos(pi/2) > 0.  So no special handling is needed for pi/2.  Even if
angles were handled in such a way that 90deg were exactly represented,
these routines would still "work", with, e.g., tan(pi/2) -> inf.

(A caution: with long doubles = a 64-bit fraction, we have cos(pi/2) <
0; and now we would need to be careful.)

As a consequence, there no need for error handling in pj_tsfn; the
HUGE_VAL return has gone and, of course, HUGE_VAL is a perfectly legal
input to tsfn's inverse, phi2, which would return -pi/2.  This "error
handling" was only needed for e = 1, a case which is filtered out
upstream.  I will note that bad argument handling is much more natural
using NAN instead of HUGE_VAL.  See issue #2376

I've renamed the error condition for non-convergence of the inverse
projection from "non-convergent inverse phi2" to "non-convergent
sinh(psi) to tan(phi)".

Now that pj_tsfn and pj_phi2 now return "better" results, there were
some malfunctions in the projections that called them, specifically
gstmerc, lcc, and tobmerc.

  * gstmerc invoked pj_tsfn(phi, sinphi, e) with a value of sinphi
    that wasn't equal to sin(phi).  Disaster followed.  I fixed this.
    I also replaced numerous occurrences of "-1.0 * x" by "-x".
    (Defining a function with arguments phi and sinphi is asking for
    trouble.)

  * lcc incorrectly thinks that the projection isn't defined for
    standard latitude = +/- 90d.  This happens to be false (it reduces
    to polar stereographic in this limit).  The check was whether
    tsfn(phi) = 0 (which only tested for the north pole not the south
    pole).  However since tsfn(pi/2) now (correctly) returns a nonzero
    result, this test fails.  I now just test for |phi| = pi/2.  This
    is clearer and catches both poles (I'm assuming that the current
    implementation will probably fail in these cases).

  * tobmerc similarly thinks that phi close to +/- pi/2 can't be
    transformed even though psi(pi/2) is only 38.  I'm disincline to
    fight this.  However I did tighten up the failure condition
    (strict equality of |phi| == pi/2).

OTHER STUFF
===========

Testing
-------

builtins.gei: I tightened up the tests for merc (and while I was about
it etmerc and tmerc) to reflect full double precision accuracy.  My
test values are generated with MPFR enabled code and so should be
accurate to all digits given.  For the record, for GRS80 I use f =
1/298.2572221008827112431628366 in these calculations.

pj_phi2_test: many of the tests were bogus testing irrelevant input
parameters, like negative values of exp(-psi), and freezing in the
arbitrary behavior of phi2.  I've reworked most for the tests to be
semi-useful.  @schwehr can you review.

Documentation
-------------

I've updated merc.rst to outline the calculation of the inverse
projection.

phi2.cpp includes detailed notes about applying Newton's method to
find tan(phi) in terms of sinh(psi).

Future work
-----------

lcc needs some tender loving care.  It can easily (and should) be
modified to allow stdlat = +/- 90 (reduces to polar stereographic),
stdlat = 0 and stdlat_1 + stdlat_2 = 0 (reduces to Mercator).  A
little more elbow grease will allow the treatment of stdlat_1 close to
stdlat_2 using divided differences.  (See my implementation of the
LambertConformalConic class in GeographicLib.)

All the places where pj_tsfn and pj_phi2 are called need to be
reworked to cut out the use of Snyder's t = exp(-psi() variable and
instead use sinh(psi).

Maybe include the machinery for series conversions between all
auxiliary latitudes as "support functions".  Then etmerc could use
this (as could mlfn for computing meridional distance).  merc could
offer the etmerc style projection via chi as an option when the
flattening is sufficiently small.

8 files changed, 168 insertions, 114 deletions

diff --git a/src/phi2.cpp b/src/phi2.cpp
index b81456b0..eb6d5c82 100644
--- a/src/phi2.cpp
+++ b/src/phi2.cpp
@@ -1,6 +1,7 @@
 /* Determine latitude angle phi-2. */
 
 #include <math.h>
+#include <limits>
 
 #include "proj.h"
 #include "proj_internal.h"
@@ -8,61 +9,125 @@
 static const double TOL = 1.0e-10;
 static const int N_ITER = 15;
 
+double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) {
+  /****************************************************************************
+  * Convert tau' = sinh(psi) = tan(chi) to tau = tan(phi).  The code is taken
+  * from GeographicLib::Math::tauf(taup, e).
+  *
+  * Here
+  *   phi = geographic latitude (radians)
+  * psi is the isometric latitude
+  *   psi = asinh(tan(phi)) - e * atanh(e * sin(phi))
+  *       = asinh(tan(chi))
+  * chi is the conformal latitude
+  *
+  * The representation of latitudes via their tangents, tan(phi) and tan(chi),
+  * maintains full *relative* accuracy close to latitude = 0 and +/- pi/2.
+  * This is sometimes important, e.g., to compute the scale of the transverse
+  * Mercator projection which involves cos(phi)/cos(chi) tan(phi)
+  *
+  * From Karney (2011), Eq. 7,
+  *
+  *   tau' = sinh(psi) = sinh(asinh(tan(phi)) - e * atanh(e * sin(phi)))
+  *        = tan(phi) * cosh(e * atanh(e * sin(phi))) -
+  *          sec(phi) * sinh(e * atanh(e * sin(phi)))
+  *        = tau * sqrt(1 + sigma^2) - sqrt(1 + tau^2) * sigma
+  * where
+  *   sigma = sinh(e * atanh( e * tau / sqrt(1 + tau^2) ))
+  *
+  * For e small,
+  *
+  *    tau' = (1 - e^2) * tau
+  *
+  * The relation tau'(tau) can therefore by reliably inverted by Newton's
+  * method with
+  *
+  *    tau = tau' / (1 - e^2)
+  *
+  * as an initial guess.  Newton's method requires dtau'/dtau.  Noting that
+  *
+  *   dsigma/dtau = e^2 * sqrt(1 + sigma^2) /
+  *                 (sqrt(1 + tau^2) * (1 + (1 - e^2) * tau^2))
+  *   d(sqrt(1 + tau^2))/dtau = tau / sqrt(1 + tau^2)
+  *
+  * we have
+  *
+  *   dtau'/dtau = (1 - e^2) * sqrt(1 + tau'^2) * sqrt(1 + tau^2) /
+  *                (1 + (1 - e^2) * tau^2)
+  *
+  * This works fine unless tau^2 and tau'^2 overflows.  This may be partially
+  * cured by writing, e.g., sqrt(1 + tau^2) as hypot(1, tau).  However, nan
+  * will still be generated with tau' = inf, since (inf - inf) will appear in
+  * the Newton iteration.
+  *
+  * If we note that for sufficiently large |tau|, i.e., |tau| >= 2/sqrt(eps),
+  * sqrt(1 + tau^2) = |tau| and
+  *
+  *   tau' = exp(- e * atanh(e)) * tau
+  *
+  * So
+  *
+  *   tau = exp(e * atanh(e)) * tau'
+  *
+  * can be returned unless |tau| >= 2/sqrt(eps); this then avoids overflow
+  * problems for large tau' and returns the correct result for tau' = +/-inf
+  * and nan.
+  *
+  * Newton's method usually take 2 iterations to converge to double precision
+  * accuracy (for WGS84 flattening).  However only 1 iteration is needed for
+  * |chi| < 3.35 deg.  In addition, only 1 iteration is needed for |chi| >
+  * 89.18 deg (tau' > 70), if tau = exp(e * atanh(e)) * tau' is used as the
+  * starting guess.
+  ****************************************************************************/
+
+  constexpr int numit = 5;
+  // min iterations = 1, max iterations = 2; mean = 1.954
+  constexpr double tol = sqrt(std::numeric_limits<double>::epsilon()) / 10;
+  constexpr double tmax = 2 / sqrt(std::numeric_limits<double>::epsilon());
+  double
+    e2m = 1 - e * e,
+    tau = fabs(taup) > 70 ? taup * exp(e * atanh(e)) : taup / e2m,
+    stol = tol * std::max(1.0, fabs(taup));
+  if (!(fabs(tau) < tmax)) return tau; // handles +/-inf and nan and e = 1
+  // if (e2m < 0) return std::numeric_limits<double>::quiet_NaN();
+  int i = numit;
+  for (; i; --i) {
+    double tau1 = sqrt(1 + tau * tau),
+      sig = sinh( e * atanh(e * tau / tau1) ),
+      taupa = sqrt(1 + sig * sig) * tau - sig * tau1,
+      dtau = (taup - taupa) * (1 + e2m * tau * tau) /
+      ( e2m * sqrt(1 + tau * tau) * sqrt(1 + taupa * taupa) );
+    tau += dtau;
+    if (!(fabs(dtau) >= stol))
+      break;
+  }
+  if (i == 0)
+    pj_ctx_set_errno(ctx, PJD_ERR_NON_CONV_SINHPSI2TANPHI);
+  return tau;
+}
+
 /*****************************************************************************/
 double pj_phi2(projCtx ctx, const double ts0, const double e) {
-/******************************************************************************
-Determine latitude angle phi-2.
-Inputs:
-  ts = exp(-psi) where psi is the isometric latitude (dimensionless)
-  e = eccentricity of the ellipsoid (dimensionless)
-Output:
-  phi = geographic latitude (radians)
-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))
-This routine inverts this relation using the iterative scheme given
-by Snyder (1987), Eqs. (7-9) - (7-11)
-*******************************************************************************/
-    const double eccnth = .5 * e;
-    double ts = ts0;
-#ifdef no_longer_used_original_convergence_on_exact_dphi
-    double Phi = M_HALFPI - 2. * atan(ts);
-#endif
-    int i = N_ITER;
-
-    for(;;) {
-        /*
-         * sin(Phi) = sin(PI/2 - 2* atan(ts))
-         *          = cos(2*atan(ts))
-         *          = 2*cos^2(atan(ts)) - 1
-         *          = 2 / (1 + ts^2) - 1
-         *          = (1 - ts^2) / (1 + ts^2)
-         */
-        const double sinPhi = (1 - ts * ts) / (1 + ts * ts);
-        const double con = e * sinPhi;
-        double old_ts = ts;
-        ts = ts0 * pow((1. - con) / (1. + con), eccnth);
-#ifdef no_longer_used_original_convergence_on_exact_dphi
-        /* The convergence criterion is nominally on exact dphi */
-        const double newPhi = M_HALFPI - 2. * atan(ts);
-        const double dphi = newPhi - Phi;
-        Phi = newPhi;
-#else
-        /* If we don't immediately update Phi from this, we can
-         * change the conversion criterion to save us computing atan() at each step.
-         * Particularly we can observe that:
-         * |atan(ts) - atan(old_ts)| <= |ts - old_ts|
-         * So if |ts - old_ts| matches our convergence criterion, we're good.
-         */
-        const double dphi = 2 * (ts - old_ts);
-#endif
-        if (fabs(dphi) > TOL && --i) {
-            continue;
-        }
-        break;
-    }
-    if (i <= 0)
-        pj_ctx_set_errno(ctx, PJD_ERR_NON_CON_INV_PHI2);
-    return M_HALFPI - 2. * atan(ts);
+  /****************************************************************************
+   * Determine latitude angle phi-2.
+   * Inputs:
+   *   ts = exp(-psi) where psi is the isometric latitude (dimensionless)
+   *   e = eccentricity of the ellipsoid (dimensionless)
+   * Output:
+   *   phi = geographic latitude (radians)
+   * 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))
+   *
+   * OLD: This routine inverts this relation using the iterative scheme given
+   * by Snyder (1987), Eqs. (7-9) - (7-11).
+   *
+   * NEW: This routine writes converts t = exp(-psi) to
+   *
+   *   tau' = sinh(psi) = (1/t - t)/2
+   *
+   * returns atan(sinpsi2tanphi(tau'))
+   ***************************************************************************/
+  return atan(pj_sinhpsi2tanphi(ctx, (1/ts0 - ts0) / 2, e));
 }
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 79b1da10..203765a3 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -626,7 +626,7 @@ struct FACTORS {
 #define PJD_ERR_INVALID_X_OR_Y          -15
 #define PJD_ERR_WRONG_FORMAT_DMS_VALUE  -16
 #define PJD_ERR_NON_CONV_INV_MERI_DIST  -17
-#define PJD_ERR_NON_CON_INV_PHI2        -18
+#define PJD_ERR_NON_CONV_SINHPSI2TANPHI -18
 #define PJD_ERR_ACOS_ASIN_ARG_TOO_LARGE -19
 #define PJD_ERR_TOLERANCE_CONDITION     -20
 #define PJD_ERR_CONIC_LAT_EQUAL         -21
@@ -846,6 +846,7 @@ double  pj_qsfn(double, double, double);
 double  pj_tsfn(double, double, double);
 double  pj_msfn(double, double, double);
 double  PROJ_DLL pj_phi2(projCtx_t *, const double, const double);
+double  pj_sinhpsi2tanphi(projCtx_t *, const double, const double);
 double  pj_qsfn_(double, PJ *);
 double *pj_authset(double);
 double  pj_authlat(double, double *);
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..b83b8072 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( abs(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( abs(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..472cb43f 100644
--- a/src/projections/merc.cpp
+++ b/src/projections/merc.cpp
@@ -10,45 +10,25 @@
 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));
+    xy.y = P->k0 * (asinh(tan(lp.phi)) - P->e * atanh(P->e * sin(lp.phi)));
     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;
 }
 
diff --git a/src/strerrno.cpp b/src/strerrno.cpp
index 5ae0d7e1..1c8673d0 100644
--- a/src/strerrno.cpp
+++ b/src/strerrno.cpp
@@ -27,7 +27,7 @@ pj_err_list[] = {
     "invalid x or y",                                                  /* -15 */
     "improperly formed DMS value",                                     /* -16 */
     "non-convergent inverse meridional dist",                          /* -17 */
-    "non-convergent inverse phi2",                                     /* -18 */
+    "non-convergent sinh(psi) to tan(phi)",                            /* -18 */
     "acos/asin: |arg| >1.+1e-14",                                      /* -19 */
     "tolerance condition error",                                       /* -20 */
     "conic lat_1 = -lat_2",                                            /* -21 */
diff --git a/src/tsfn.cpp b/src/tsfn.cpp
index 32da09f2..fe8f29ed 100644
--- a/src/tsfn.cpp
+++ b/src/tsfn.cpp
@@ -4,14 +4,25 @@
 #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)
+   * 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))
+   ***************************************************************************/
 
-    /* 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 );
 }