diff options
Diffstat (limited to 'src/phi2.cpp')
| -rw-r--r-- | src/phi2.cpp | 45 |
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); } |