diff options
| -rw-r--r-- | src/phi2.cpp | 23 | ||||
| -rw-r--r-- | src/projections/merc.cpp | 3 |
2 files changed, 12 insertions, 14 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; diff --git a/src/projections/merc.cpp b/src/projections/merc.cpp index 39685478..3a0ed7b4 100644 --- a/src/projections/merc.cpp +++ b/src/projections/merc.cpp @@ -15,7 +15,8 @@ static PJ_XY merc_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward xy.x = P->k0 * lp.lam; // Instead of calling tan and sin, call sin and cos which the compiler // optimizes to a single call to sincos. - double sphi = sin(lp.phi), cphi = cos(lp.phi); + double sphi = sin(lp.phi); + double cphi = cos(lp.phi); xy.y = P->k0 * (asinh(sphi/cphi) - P->e * atanh(P->e * sphi)); return xy; } |