oxiproj-core 0.1.2

//! Meridional-arc length functions ported from PROJ `src/mlfn.cpp` (Karney/auxiliary-latitude based).
//!
//! These compute the meridional (meridian-arc) distance for an ellipsoid and
//! its inverse using a 6th-order expansion in the third flattening `n`. This
//! gives full double-precision accuracy for `|f| <= 1/150`. The heavy lifting
//! (Fourier-coefficient generation and Clenshaw evaluation) lives in the
//! auxiliary-latitude module [`crate::math::aux`]; this module wires those
//! pieces together exactly as `src/mlfn.cpp` does.

use crate::error::ProjResult;

/// Ported from src/mlfn.cpp (the `en` coefficient array returned by `pj_enfn`).
///
/// `data[0]` is the rectifying radius; `data[1..7]` are the geographic→rectifying
/// Fourier coefficients; `data[7..13]` are the rectifying→geographic coefficients.
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct EnFn {
    /// The 13-element coefficient table. In PROJ this is a malloc'd
    /// `double[2*Lmax + 1]` (with `Lmax == 6`), so `2*6 + 1 == 13` doubles.
    pub data: [f64; 13],
}

/// Ported from src/mlfn.cpp (`pj_enfn`).
///
/// Builds the meridional-distance coefficient table for the given third
/// flattening `n`. `data[0]` receives the rectifying radius; `data[1..7]`
/// receive the geographic→rectifying coefficients; `data[7..13]` receive the
/// rectifying→geographic coefficients. The two `&mut data[..]` borrows are
/// sequential (non-overlapping), so there is no borrow conflict, and the
/// constant ranges are in bounds for a length-13 array.
pub fn pj_enfn(n: f64) -> ProjResult<EnFn> {
    use crate::math::aux::{pj_auxlat_coeffs, pj_rectifying_radius, AuxLat};
    let mut data = [0.0_f64; 13];
    data[0] = pj_rectifying_radius(n);
    pj_auxlat_coeffs(n, AuxLat::Geographic, AuxLat::Rectifying, &mut data[1..7])?;
    pj_auxlat_coeffs(n, AuxLat::Rectifying, AuxLat::Geographic, &mut data[7..13])?;
    Ok(EnFn { data })
}

/// Ported from src/mlfn.cpp (`pj_mlfn`).
///
/// Computes the meridional distance for geographic latitude `phi` (given its
/// sine `sphi` and cosine `cphi`) using the coefficient table `en`. Mirrors
/// `en[0] * pj_auxlat_convert(phi, sphi, cphi, en + 1)`.
pub fn pj_mlfn(phi: f64, sphi: f64, cphi: f64, en: &EnFn) -> f64 {
    use crate::math::aux::pj_auxlat_convert_sc;
    en.data[0] * pj_auxlat_convert_sc(phi, sphi, cphi, &en.data[1..7], 6)
}

/// Ported from src/mlfn.cpp (`pj_inv_mlfn`).
///
/// Inverts [`pj_mlfn`]: given a meridional distance `mu`, recovers the
/// geographic latitude using the rectifying→geographic coefficients in `en`.
/// Mirrors `pj_auxlat_convert(mu / en[0], en + 1 + Lmax)`.
pub fn pj_inv_mlfn(mu: f64, en: &EnFn) -> f64 {
    use crate::math::aux::pj_auxlat_convert;
    pj_auxlat_convert(mu / en.data[0], &en.data[7..13], 6)
}

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

    /// WGS84 third flattening `n = f / (2 - f)`.
    fn wgs84_n() -> f64 {
        let f = 1.0 / 298.257223563_f64;
        f / (2.0 - f)
    }

    #[test]
    fn round_trip_sample_latitudes() {
        let n = wgs84_n();
        let en = pj_enfn(n).expect("enfn");
        for &phi in &[-1.4, -0.7, -0.1, 0.0, 0.1, 0.7, 1.4] {
            let m = pj_mlfn(phi, phi.sin(), phi.cos(), &en);
            let back = pj_inv_mlfn(m, &en);
            assert!((back - phi).abs() < 1e-12, "phi={} back={}", phi, back);
        }
    }

    #[test]
    fn round_trip_dense_sweep() {
        let n = wgs84_n();
        let en = pj_enfn(n).expect("enfn");
        let mut phi = -1.45_f64;
        while phi <= 1.45 {
            let m = pj_mlfn(phi, phi.sin(), phi.cos(), &en);
            let back = pj_inv_mlfn(m, &en);
            assert!((back - phi).abs() < 1e-11, "phi={} back={}", phi, back);
            phi += 0.02;
        }
    }

    #[test]
    fn mlfn_at_equator_is_zero() {
        let n = wgs84_n();
        let en = pj_enfn(n).expect("enfn");
        let m = pj_mlfn(0.0, 0.0, 1.0, &en);
        assert!(m.abs() < 1e-15, "m={}", m);
    }

    #[test]
    fn rectifying_radius_in_expected_range() {
        let n = wgs84_n();
        let en = pj_enfn(n).expect("enfn");
        // For WGS84's small positive n, the rectifying radius is slightly < 1.
        assert!(en.data[0] > 0.99, "radius={}", en.data[0]);
        assert!(en.data[0] < 1.0, "radius={}", en.data[0]);
    }

    #[test]
    fn enfn_layout_radius_then_fwd_then_inv() {
        use crate::math::aux::{pj_auxlat_coeffs, pj_rectifying_radius, AuxLat};
        let n = wgs84_n();
        let en = pj_enfn(n).expect("enfn");
        // data[0] must equal the rectifying radius.
        assert_eq!(en.data[0], pj_rectifying_radius(n));
        // data[1..7] must equal the geographic->rectifying coefficients.
        let mut fwd = [0.0_f64; 6];
        pj_auxlat_coeffs(n, AuxLat::Geographic, AuxLat::Rectifying, &mut fwd).expect("fwd coeffs");
        assert_eq!(&en.data[1..7], &fwd[..]);
        // data[7..13] must equal the rectifying->geographic coefficients.
        let mut inv = [0.0_f64; 6];
        pj_auxlat_coeffs(n, AuxLat::Rectifying, AuxLat::Geographic, &mut inv).expect("inv coeffs");
        assert_eq!(&en.data[7..13], &inv[..]);
    }
}