mgrs-convert 0.1.0

//! WGS84 constants and the conformal-latitude helpers shared by both projections.
//!
//! The two geodetic parameters below are the only ones in the crate, each pinned
//! by a unit test. The tau functions are the conformal mapping from Karney,
//! "Transverse Mercator with an accuracy of a few nanometers" (2011), section 3.
//
// This module performs bounded floating-point arithmetic only, which cannot
// side-effect, so `clippy::arithmetic_side_effects` is not denied here. See the
// errors section of the behavioral specification.

/// WGS84 semi-major axis in metres, exact.
pub(crate) const A: f64 = 6_378_137.0;

/// WGS84 flattening, exact.
pub(crate) const F: f64 = 1.0 / 298.257_223_563;

/// First eccentricity squared, e^2 = f(2 - f).
pub(crate) const E2: f64 = F * (2.0 - F);

/// First eccentricity.
pub(crate) fn es() -> f64 {
    E2.sqrt()
}

/// tau' = tan of the conformal latitude as a function of tau = tan(phi).
pub(crate) fn taupf(tau: f64) -> f64 {
    let tau1 = (1.0 + tau * tau).sqrt();
    let sig = (es() * (es() * tau / tau1).atanh()).sinh();
    (1.0 + sig * sig).sqrt() * tau - sig * tau1
}

/// Inverse of [`taupf`] by Newton's method, following the reference paper.
pub(crate) fn tauf(taup: f64) -> f64 {
    let e2m = 1.0 - E2;
    let mut tau = taup / e2m;
    let stol = 0.1 * f64::EPSILON.sqrt() * taup.abs().max(1.0);
    // Five iterations converge to full double precision for every finite input;
    // the loop bound keeps the function total.
    for _ in 0..5 {
        let taupa = taupf(tau);
        let dtau = (taup - taupa) * (1.0 + e2m * tau * tau)
            / (e2m * (1.0 + tau * tau).sqrt() * (1.0 + taupa * taupa).sqrt());
        tau += dtau;
        if dtau.abs() < stol {
            break;
        }
    }
    tau
}