aboutsummaryrefslogtreecommitdiff
path: root/src/phi2.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/phi2.cpp')
-rw-r--r--src/phi2.cpp175
1 files changed, 120 insertions, 55 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));
}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 21:28:16 +0000