pj_phi2(): speed-up computation (and thus inverse ellipsoidal Mercator and LCC)

This does not change the numeric values returned by the method,
as far as I could see on a few samplings.

The tricks used save a call to sin() and atan() at each iteration.

This directly affects speed of inverse Mercator and LCC (among others),
in their ellipsoidal formulation.

Timings on inverse Mercator show a 31% speed-up at mid-latitudes
where pj_phi2() needs 5 iterations, and 24% at latitudes close to 0 or
90deg where it needs one iteration.

1 files changed, 33 insertions, 12 deletions

diff --git a/src/phi2.cpp b/src/phi2.cpp
index afabd06e..b81456b0 100644
--- a/src/phi2.cpp
+++ b/src/phi2.cpp
@@ -9,7 +9,7 @@ static const double TOL = 1.0e-10;
 static const int N_ITER = 15;
 
 /*****************************************************************************/
-double pj_phi2(projCtx ctx, double ts, double e) {
+double pj_phi2(projCtx ctx, const double ts0, const double e) {
 /******************************************************************************
 Determine latitude angle phi-2.
 Inputs:
@@ -24,24 +24,45 @@ Here isometric latitude is defined by
 This routine inverts this relation using the iterative scheme given
 by Snyder (1987), Eqs. (7-9) - (7-11)
 *******************************************************************************/
-    double eccnth = .5 * e;
+    const double eccnth = .5 * e;
+    double ts = ts0;
+#ifdef no_longer_used_original_convergence_on_exact_dphi
     double Phi = M_HALFPI - 2. * atan(ts);
-    double con;
+#endif
     int i = N_ITER;
 
     for(;;) {
-        double dphi;
-        con = e * sin(Phi);
-        dphi = M_HALFPI - 2. * atan(ts * pow((1. - con) /
-           (1. + con), eccnth)) - Phi;
-
-        Phi += dphi;
-
-        if (fabs(dphi) > TOL && --i)
+        /*
+         * 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 Phi;
+    return M_HALFPI - 2. * atan(ts);
 }