oxiproj-core 0.1.2

//! `pj_phi2` / `pj_sinhpsi2tanphi` ported from PROJ `src/phi2.cpp` (GeographicLib `tauf`).
//!
//! These helpers invert the relationship between the isometric latitude `psi`
//! (equivalently `ts = exp(-psi)` or `taup = sinh(psi) = tan(chi)`, where `chi`
//! is the conformal latitude) and the geographic latitude `phi`. The inversion
//! uses Newton's method on `tau = tan(phi)`, following
//! `GeographicLib::Math::tauf(taup, e)` as ported into PROJ's `pj_sinhpsi2tanphi`.
//!
//! Faithful note: PROJ records a context errno (`PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE`)
//! when the Newton iteration fails to converge within `numit` steps; this port
//! has no projection context at this layer, so on non-convergence it simply
//! returns the last `tau` estimate (mirroring PROJ's numerical behaviour without
//! the side-effecting error report).

/// Ported from src/phi2.cpp (`pj_sinhpsi2tanphi`).
///
/// Given `taup = sinh(psi)` returns `tan(phi)` via Newton's method (GeographicLib).
///
/// Faithful note: PROJ records a context errno when the iteration fails to
/// converge within `numit` steps; this port has no context, so on
/// non-convergence it simply returns the last `tau` estimate.
pub fn pj_sinhpsi2tanphi(taup: f64, e: f64) -> f64 {
    const NUMIT: i32 = 5;
    let rooteps = f64::EPSILON.sqrt();
    let tol = rooteps / 10.0;
    let tmax = 2.0 / rooteps;
    let e2m = 1.0 - e * e;
    let stol = tol * taup.abs().max(1.0);
    // Initial guess. 70 corresponds to chi = 89.18 deg (see PROJ commentary).
    let mut tau = if taup.abs() > 70.0 {
        taup * (e * e.atanh()).exp()
    } else {
        taup / e2m
    };
    // Handles +/-inf and nan and e = 1 (matches `!(fabs(tau) < tmax)` in C).
    #[allow(clippy::neg_cmp_op_on_partial_ord)]
    if !(tau.abs() < tmax) {
        return tau;
    }
    // C loop: `for (i = numit; i; --i)`. The body runs while `i != 0`, the
    // convergence break is tested before the decrement, so the decrement here
    // happens AFTER the convergence check, giving at most NUMIT iterations.
    let mut i = NUMIT;
    while i > 0 {
        let tau1 = (1.0 + tau * tau).sqrt();
        let sig = (e * (e * tau / tau1).atanh()).sinh();
        let taupa = (1.0 + sig * sig).sqrt() * tau - sig * tau1;
        let dtau = ((taup - taupa) * (1.0 + e2m * tau * tau))
            / (e2m * tau1 * (1.0 + taupa * taupa).sqrt());
        tau += dtau;
        // Backwards test to allow nans to succeed (mirrors `!(fabs(dtau) >= stol)`).
        #[allow(clippy::neg_cmp_op_on_partial_ord)]
        if !(dtau.abs() >= stol) {
            break;
        }
        i -= 1;
    }
    tau
}

/// Ported from src/phi2.cpp (`pj_phi2`).
///
/// Converts `ts0 = exp(-psi)` (the variable `ts` of Snyder (1987), Eq. (7-10))
/// to `tau' = tan(chi) = sinh(psi) = (1/ts0 - ts0) / 2`, then returns
/// `atan(pj_sinhpsi2tanphi(tau'))`, i.e. the geographic latitude `phi`.
pub fn pj_phi2(ts0: f64, e: f64) -> f64 {
    pj_sinhpsi2tanphi((1.0 / ts0 - ts0) / 2.0, e).atan()
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::math::tsfn::pj_tsfn;

    /// WGS84 eccentricity, derived from the first eccentricity squared
    /// `es = 0.0066943799901413165`.
    const WGS84_E: f64 = 0.0818191908426215;

    #[test]
    fn wgs84_eccentricity_matches_es_sqrt() {
        let es = 0.0066943799901413165_f64;
        let e = es.sqrt();
        assert!((e - WGS84_E).abs() < 1e-15, "e = {e}");
    }

    #[test]
    fn sinhpsi2tanphi_zero_maps_to_zero() {
        let e = WGS84_E;
        let tau = pj_sinhpsi2tanphi(0.0, e);
        assert!(tau.abs() < 1e-15, "expected ~0, got {tau}");
    }

    #[test]
    fn phi2_at_equator_is_zero() {
        // At phi = 0, ts = pj_tsfn(0, 0, e) = 1, and pj_phi2(1, e) = 0.
        let e = WGS84_E;
        let ts = pj_tsfn(0.0, 0.0_f64.sin(), e);
        assert!((ts - 1.0).abs() < 1e-15, "expected ts ~1, got {ts}");
        let phi = pj_phi2(ts, e);
        assert!(phi.abs() < 1e-15, "expected phi ~0, got {phi}");
    }

    #[test]
    fn round_trip_selected_latitudes() {
        let e = WGS84_E;
        let phis = [-1.4_f64, -1.0, -0.5, -0.1, 0.0, 0.1, 0.5, 1.0, 1.4];
        for &phi in &phis {
            let ts = pj_tsfn(phi, phi.sin(), e);
            let recovered = pj_phi2(ts, e);
            assert!(
                (recovered - phi).abs() < 1e-12,
                "round trip failed: phi = {phi}, recovered = {recovered}, ts = {ts}"
            );
        }
    }

    #[test]
    fn round_trip_dense_sweep() {
        // proptest-like coverage without adding a dependency: sweep phi across a
        // dense grid spanning nearly the full open interval (-pi/2, pi/2).
        let e = WGS84_E;
        let mut steps = 0_u32;
        let mut phi = -1.45_f64;
        while phi <= 1.45 {
            let ts = pj_tsfn(phi, phi.sin(), e);
            let recovered = pj_phi2(ts, e);
            assert!(
                (recovered - phi).abs() < 1e-11,
                "dense round trip failed: phi = {phi}, recovered = {recovered}, ts = {ts}"
            );
            steps += 1;
            phi += 0.01;
        }
        // Sanity: ensure the loop actually exercised a wide range of values.
        assert!(
            steps > 250,
            "expected a dense sweep, only ran {steps} steps"
        );
    }
}