Function sinhpsi2tanphi 

pub fn sinhpsi2tanphi(taup: f64, e: f64) -> f64

Convert tau’ = sinh(psi) = tan(chi) to tau = tan(phi). The code is taken from GeographicLib::Math::tauf(taup, e).

Here phi = geographic latitude (radians) psi is the isometric latitude psi = asinh(tan(phi)) - e * atanh(e * sin(phi)) = asinh(tan(chi)) chi is the conformal latitude

The representation of latitudes via their tangents, tan(phi) and tan(chi), maintains full relative accuracy close to latitude = 0 and +/- pi/2. This is sometimes important, e.g., to compute the scale of the transverse Mercator projection which involves cos(phi)/cos(chi) tan(phi)

From Karney (2011), Eq. 7,

tau’ = sinh(psi) = sinh(asinh(tan(phi)) - e * atanh(e * sin(phi))) = tan(phi) * cosh(e * atanh(e * sin(phi))) - sec(phi) * sinh(e * atanh(e * sin(phi))) = tau * sqrt(1 + sigma^2) - sqrt(1 + tau^2) * sigma where sigma = sinh(e * atanh( e * tau / sqrt(1 + tau^2) ))

For e small,

tau’ = (1 - e^2) * tau

The relation tau’(tau) can therefore by reliably inverted by Newton’s method with

tau = tau’ / (1 - e^2)

as an initial guess. Newton’s method requires dtau’/dtau. Noting that

dsigma/dtau = e^2 * sqrt(1 + sigma^2) / (sqrt(1 + tau^2) * (1 + (1 - e^2) * tau^2)) d(sqrt(1 + tau^2))/dtau = tau / sqrt(1 + tau^2)

we have

dtau’/dtau = (1 - e^2) * sqrt(1 + tau’^2) * sqrt(1 + tau^2) / (1 + (1 - e^2) * tau^2)

This works fine unless tau^2 and tau’^2 overflows. This may be partially cured by writing, e.g., sqrt(1 + tau^2) as hypot(1, tau). However, nan will still be generated with tau’ = inf, since (inf - inf) will appear in the Newton iteration.

If we note that for sufficiently large |tau|, i.e., |tau| >= 2/sqrt(eps), sqrt(1 + tau^2) = |tau| and

tau’ = exp(- e * atanh(e)) * tau

So

tau = exp(e * atanh(e)) * tau’

can be returned unless |tau| >= 2/sqrt(eps); this then avoids overflow problems for large tau’ and returns the correct result for tau’ = +/-inf and nan.

Newton’s method usually take 2 iterations to converge to double precision accuracy (for WGS84 flattening). However only 1 iteration is needed for |chi| < 3.35 deg. In addition, only 1 iteration is needed for |chi| > 89.18 deg (tau’ > 70), if tau = exp(e * atanh(e)) * tau’ is used as the starting guess.

For small flattening, |f| <= 1/50, the series expansion in n can be used:

Assuming n = e^2 / (1 + sqrt(1 - e^2))^2 is passed as an argument

double F[int(AuxLat::AUXORDER)]; auxlat_coeffs(n, AuxLat::CONFORMAL, AuxLat::GEOGRAPHIC, F); double sphi, cphi; // schi cchi auxlat_convert(taup/hypot(1.0, taup), 1/hypot(1.0, taup), sphi, cphi, F); return sphi/cphi;