diff options
Diffstat (limited to 'src/phi2.cpp')
| -rw-r--r-- | src/phi2.cpp | 23 |
1 files changed, 10 insertions, 13 deletions
diff --git a/src/phi2.cpp b/src/phi2.cpp index 0fdca47c..1c48d67f 100644 --- a/src/phi2.cpp +++ b/src/phi2.cpp @@ -80,13 +80,11 @@ double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) { constexpr int numit = 5; // min iterations = 1, max iterations = 2; mean = 1.954 - static const double - rooteps = sqrt(std::numeric_limits<double>::epsilon()), - tol = rooteps / 10, // the convergence criterion for Newton's method - tmax = 2 / rooteps; // the large arg limit is exact for tau > tmax - const double - e2m = 1 - e * e, - stol = tol * std::max(1.0, fabs(taup)); + static const double rooteps = sqrt(std::numeric_limits<double>::epsilon()); + static const double tol = rooteps / 10; // the criterion for Newton's method + static const double tmax = 2 / rooteps; // threshold for large arg limit exact + const double e2m = 1 - e * e; + const double stol = tol * std::max(1.0, fabs(taup)); // The initial guess. 70 corresponds to chi = 89.18 deg (see above) double tau = fabs(taup) > 70 ? taup * exp(e * atanh(e)) : taup / e2m; if (!(fabs(tau) < tmax)) // handles +/-inf and nan and e = 1 @@ -95,12 +93,11 @@ double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) { // 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 * tau1 * sqrt(1 + taupa * taupa)) ); + double tau1 = sqrt(1 + tau * tau); + double sig = sinh( e * atanh(e * tau / tau1) ); + double taupa = sqrt(1 + sig * sig) * tau - sig * tau1; + double dtau = ( (taup - taupa) * (1 + e2m * (tau * tau)) / + (e2m * tau1 * sqrt(1 + taupa * taupa)) ); tau += dtau; if (!(fabs(dtau) >= stol)) // backwards test to allow nans to succeed. break; |