| 8 | #include "proj_internal.h" |
| 9 | |
| 10 | double pj_sinhpsi2tanphi(PJ_CONTEXT *ctx, const double taup, const double e) { |
| 11 | /*************************************************************************** |
| 12 | * Convert tau' = sinh(psi) = tan(chi) to tau = tan(phi). The code is taken |
| 13 | * from GeographicLib::Math::tauf(taup, e). |
| 14 | * |
| 15 | * Here |
| 16 | * phi = geographic latitude (radians) |
| 17 | * psi is the isometric latitude |
| 18 | * psi = asinh(tan(phi)) - e * atanh(e * sin(phi)) |
| 19 | * = asinh(tan(chi)) |
| 20 | * chi is the conformal latitude |
| 21 | * |
| 22 | * The representation of latitudes via their tangents, tan(phi) and |
| 23 | * tan(chi), maintains full *relative* accuracy close to latitude = 0 and |
| 24 | * +/- pi/2. This is sometimes important, e.g., to compute the scale of the |
| 25 | * transverse Mercator projection which involves cos(phi)/cos(chi) tan(phi) |
| 26 | * |
| 27 | * From Karney (2011), Eq. 7, |
| 28 | * |
| 29 | * tau' = sinh(psi) = sinh(asinh(tan(phi)) - e * atanh(e * sin(phi))) |
| 30 | * = tan(phi) * cosh(e * atanh(e * sin(phi))) - |
| 31 | * sec(phi) * sinh(e * atanh(e * sin(phi))) |
| 32 | * = tau * sqrt(1 + sigma^2) - sqrt(1 + tau^2) * sigma |
| 33 | * where |
| 34 | * sigma = sinh(e * atanh( e * tau / sqrt(1 + tau^2) )) |
| 35 | * |
| 36 | * For e small, |
| 37 | * |
| 38 | * tau' = (1 - e^2) * tau |
| 39 | * |
| 40 | * The relation tau'(tau) can therefore by reliably inverted by Newton's |
| 41 | * method with |
| 42 | * |
| 43 | * tau = tau' / (1 - e^2) |
| 44 | * |
| 45 | * as an initial guess. Newton's method requires dtau'/dtau. Noting that |
| 46 | * |
| 47 | * dsigma/dtau = e^2 * sqrt(1 + sigma^2) / |
| 48 | * (sqrt(1 + tau^2) * (1 + (1 - e^2) * tau^2)) |
| 49 | * d(sqrt(1 + tau^2))/dtau = tau / sqrt(1 + tau^2) |
| 50 | * |
| 51 | * we have |
| 52 | * |
| 53 | * dtau'/dtau = (1 - e^2) * sqrt(1 + tau'^2) * sqrt(1 + tau^2) / |
| 54 | * (1 + (1 - e^2) * tau^2) |
| 55 | * |
| 56 | * This works fine unless tau^2 and tau'^2 overflows. This may be partially |
| 57 | * cured by writing, e.g., sqrt(1 + tau^2) as hypot(1, tau). However, nan |
| 58 | * will still be generated with tau' = inf, since (inf - inf) will appear in |
| 59 | * the Newton iteration. |
| 60 | * |
| 61 | * If we note that for sufficiently large |tau|, i.e., |tau| >= 2/sqrt(eps), |
| 62 | * sqrt(1 + tau^2) = |tau| and |
| 63 | * |
| 64 | * tau' = exp(- e * atanh(e)) * tau |
| 65 | * |
| 66 | * So |
| 67 | * |
no test coverage detected