diff options
Diffstat (limited to 'src/tsfn.cpp')
| -rw-r--r-- | src/tsfn.cpp | 29 |
1 files changed, 20 insertions, 9 deletions
diff --git a/src/tsfn.cpp b/src/tsfn.cpp index 32da09f2..fe8f29ed 100644 --- a/src/tsfn.cpp +++ b/src/tsfn.cpp @@ -4,14 +4,25 @@ #include "proj_internal.h" double pj_tsfn(double phi, double sinphi, double e) { - double denominator; - sinphi *= e; + /**************************************************************************** + * Determine function ts(phi) defined in Snyder (1987), Eq. (7-10) + * Inputs: + * phi = geographic latitude (radians) + * e = eccentricity of the ellipsoid (dimensionless) + * Output: + * ts = exp(-psi) where psi is the isometric latitude (dimensionless) + * 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)) + ***************************************************************************/ - /* avoid zero division, fail gracefully */ - denominator = 1.0 + sinphi; - if (denominator == 0.0) - return HUGE_VAL; - - return (tan (.5 * (M_HALFPI - phi)) / - pow((1. - sinphi) / (denominator), .5 * e)); + double cosphi = cos(phi); + // exp(-asinh(tan(phi))) = 1 / (tan(phi) + sec(phi)) + // = cos(phi) / (1 + sin(phi)) good for phi > 0 + // = (1 - sin(phi)) / cos(phi) good for phi < 0 + return exp(e * atanh(e * sinphi)) * + ( sinphi > 0 ? + cosphi / (1 + sinphi) : + (1 - sinphi) / cosphi ); } |