oxiproj-core 0.1.2

//! Meridional distance series ported from PROJ `src/proj_mdist.cpp` (Gerald Evenden).
//!
//! Computes the distance from the equator along the meridian to a latitude
//! `phi` on the unit ellipsoid, plus its inverse, to a precision commensurate
//! with double precision. This is the classic Gerald Evenden series expansion;
//! see `src/proj_mdist.cpp` in PROJ 9.8.0.

#[cfg(feature = "no_std")]
use alloc::vec;
#[cfg(feature = "no_std")]
use alloc::vec::Vec;

/// Ported from src/proj_mdist.cpp (the `MDIST` struct).
#[derive(Debug, Clone, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct MDist {
    /// number of series terms (C `nb`).
    nb: usize,
    /// eccentricity squared (C `es`).
    es: f64,
    /// the C field `b->E` (capital E) — leading coefficient.
    big_e: f64,
    /// series coefficients (C `b->b[...]`), length `nb + 1`.
    b: Vec<f64>,
}

/// Ported from src/proj_mdist.cpp (`proj_mdist_ini`).
///
/// Builds the meridional-distance series coefficients for eccentricity-squared
/// `es`. Mirrors the C initialiser exactly, including the iteration bookkeeping:
/// the term loop runs `i = 1 .. MAX_ITER`, breaking early once the accumulated
/// `Es` stops changing. As in the C source, `E[0]` is initialised to `1.0`
/// (`double E[MAX_ITER] = {1.};`), the update `denf *= ++denfi;` is a
/// pre-increment, `nb = i - 1`, and the coefficient vector has length `i`
/// (indices `0..=nb`).
pub fn proj_mdist_ini(es: f64) -> MDist {
    const MAX_ITER: usize = 20;
    // C: `double E[MAX_ITER] = {1.};` — E[0] == 1.0, the rest 0.0.
    let mut e_series = [0.0_f64; MAX_ITER];
    e_series[0] = 1.0;
    let mut ens = es;
    let mut numf = 1.0_f64;
    let mut twon1 = 1.0_f64;
    let mut denfi = 1.0_f64;
    let mut denf = 1.0_f64;
    let mut twon = 4.0_f64;
    let mut es_acc = 1.0_f64; // C `Es`
    let mut el = 1.0_f64; // C `El`
                          // `last_i` tracks the value of C's `i` once the loop ends: the break value
                          // if it breaks, otherwise MAX_ITER after natural completion.
    let mut last_i = MAX_ITER;
    let mut i = 1usize;
    while i < MAX_ITER {
        numf *= twon1 * twon1;
        let den = twon * denf * denf * twon1;
        let t = numf / den;
        let term = t * ens;
        e_series[i] = term;
        es_acc -= term;
        ens *= es;
        twon *= 4.0;
        denfi += 1.0; // C `++denfi`
        denf *= denfi; // C `denf *= ++denfi`
        twon1 += 2.0;
        if es_acc == el {
            // C: `break;` — `i` keeps its current value (no `++i`).
            last_i = i;
            break;
        }
        el = es_acc;
        i += 1;
    }
    // If the loop completed without breaking, C's `i == MAX_ITER`.
    if i >= MAX_ITER {
        last_i = MAX_ITER;
    }
    let nb = last_i - 1;
    let big_e = es_acc;
    // C allocates `i` usable doubles for `b->b`; indices 0..=last_i-1 are filled.
    let mut b = vec![0.0_f64; last_i];
    b[0] = 1.0 - es_acc;
    let mut es2 = 1.0 - es_acc; // C reuses `Es` here
    let mut numf2 = 1.0_f64;
    let mut denf2 = 1.0_f64;
    let mut numfi = 2.0_f64;
    let mut denfi2 = 3.0_f64;
    let mut j = 1usize;
    while j < last_i {
        es2 -= e_series[j];
        numf2 *= numfi;
        denf2 *= denfi2;
        b[j] = es2 * numf2 / denf2;
        numfi += 2.0;
        denfi2 += 2.0;
        j += 1;
    }
    MDist { nb, es, big_e, b }
}

/// Ported from src/proj_mdist.cpp (`proj_mdist`).
///
/// Evaluates the meridional distance to latitude `phi` (with `sphi = sin(phi)`,
/// `cphi = cos(phi)`) using the precomputed series in `m`.
pub fn proj_mdist(phi: f64, sphi: f64, cphi: f64, m: &MDist) -> f64 {
    let sc = sphi * cphi;
    let sphi2 = sphi * sphi;
    let d = phi * m.big_e - m.es * sc / (1.0 - m.es * sphi2).sqrt();
    // C: `sum = b->b[i = b->nb]; while (i) sum = b->b[--i] + sphi2 * sum;`
    let mut sum = match m.b.get(m.nb) {
        Some(v) => *v,
        None => 0.0,
    };
    let mut i = m.nb;
    while i > 0 {
        i -= 1;
        let bi = match m.b.get(i) {
            Some(v) => *v,
            None => 0.0,
        };
        sum = bi + sphi2 * sum;
    }
    d + sc * sum
}

/// Ported from src/proj_mdist.cpp (`proj_inv_mdist`).
///
/// Inverts [`proj_mdist`] via Newton's method to recover `phi` from a
/// meridional distance `dist`.
///
/// Faithful note: PROJ records a context errno
/// (`PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN`) on non-convergence;
/// this port has no context and simply returns the last `phi` estimate after
/// `MAX_ITER` iterations (mirroring PROJ's numerical behaviour without the
/// side-effecting error report).
pub fn proj_inv_mdist(dist: f64, m: &MDist) -> f64 {
    const TOL: f64 = 1e-14;
    const MAX_ITER: usize = 20;
    let k = 1.0 / (1.0 - m.es);
    let mut phi = dist;
    let mut iter = MAX_ITER;
    while iter > 0 {
        let s = phi.sin();
        let t = 1.0 - m.es * s * s;
        let delta = (proj_mdist(phi, s, phi.cos(), m) - dist) * (t * t.sqrt()) * k;
        phi -= delta;
        if delta.abs() < TOL {
            return phi;
        }
        iter -= 1;
    }
    phi
}

#[cfg(test)]
mod tests {
    use super::*;

    /// WGS84 eccentricity squared.
    const WGS84_ES: f64 = 0.0066943799901413165;

    #[test]
    fn ini_produces_consistent_series_lengths() {
        let m = proj_mdist_ini(WGS84_ES);
        assert!(m.nb >= 1, "nb should be at least 1, got {}", m.nb);
        assert_eq!(
            m.b.len(),
            m.nb + 1,
            "b vector length should be nb + 1 (len={}, nb={})",
            m.b.len(),
            m.nb
        );
        assert_eq!(m.es, WGS84_ES);
        assert!(m.big_e.is_finite());
    }

    #[test]
    fn mdist_at_equator_is_zero() {
        let m = proj_mdist_ini(WGS84_ES);
        let d = proj_mdist(0.0, 0.0, 1.0, &m);
        assert!(
            d.abs() < 1e-15,
            "distance at equator should be ~0, got {}",
            d
        );
    }

    #[test]
    fn round_trip_sampled_latitudes() {
        let m = proj_mdist_ini(WGS84_ES);
        let phis = [-1.4_f64, -0.7, -0.1, 0.0, 0.1, 0.7, 1.4];
        for &phi in &phis {
            let d = proj_mdist(phi, phi.sin(), phi.cos(), &m);
            let back = proj_inv_mdist(d, &m);
            assert!((back - phi).abs() < 1e-11, "phi={} back={}", phi, back);
        }
    }

    #[test]
    fn round_trip_dense_sweep() {
        let m = proj_mdist_ini(WGS84_ES);
        let mut step = 0;
        loop {
            let phi = -1.45_f64 + 0.05 * (step as f64);
            if phi > 1.45 + 1e-9 {
                break;
            }
            let d = proj_mdist(phi, phi.sin(), phi.cos(), &m);
            let back = proj_inv_mdist(d, &m);
            assert!((back - phi).abs() < 1e-11, "phi={} back={}", phi, back);
            step += 1;
        }
    }

    #[test]
    fn mdist_is_monotonic_increasing() {
        let m = proj_mdist_ini(WGS84_ES);
        let lo = proj_mdist(0.1, 0.1_f64.sin(), 0.1_f64.cos(), &m);
        let hi = proj_mdist(0.7, 0.7_f64.sin(), 0.7_f64.cos(), &m);
        assert!(
            hi > lo,
            "distance should grow with latitude: lo={} hi={}",
            lo,
            hi
        );
    }

    #[test]
    fn mdist_is_odd_in_phi() {
        let m = proj_mdist_ini(WGS84_ES);
        let phi = 0.6_f64;
        let pos = proj_mdist(phi, phi.sin(), phi.cos(), &m);
        let neg = proj_mdist(-phi, (-phi).sin(), (-phi).cos(), &m);
        assert!(
            (pos + neg).abs() < 1e-13,
            "mdist should be odd: pos={} neg={}",
            pos,
            neg
        );
    }
}