oxiproj-core 0.1.2

Foundation types for OxiProj: coordinates, errors, ellipsoids, datums, and units.
Documentation
//! `pj_msfn` ported from PROJ `src/msfn.cpp`.

use crate::scalar::Scalar;

/// Ported from src/msfn.cpp (`pj_msfn`).
///
/// Computes `m = cos(phi) / sqrt(1 - es*sin^2(phi))`.
pub fn pj_msfn(sinphi: f64, cosphi: f64, es: f64) -> f64 {
    cosphi / (1.0 - es * sinphi * sinphi).sqrt()
}

/// Generic version of [`pj_msfn`] for automatic-differentiation support.
///
/// Ellipsoid constants (`es`) stay as `f64`; only coordinate values become
/// generic over `S: Scalar`. At `S = f64` this is identical to [`pj_msfn`].
pub fn pj_msfn_g<S: Scalar>(sinphi: S, cosphi: S, es: f64) -> S {
    let es_s = S::from_f64(es);
    cosphi / (S::one() - es_s * sinphi * sinphi).sqrt()
}

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

    #[test]
    fn msfn_sphere_equals_cosphi() {
        let phi = 0.5_f64;
        let m = pj_msfn(phi.sin(), phi.cos(), 0.0);
        assert!((m - phi.cos()).abs() < 1e-12);
    }

    #[test]
    fn msfn_ellipsoid_is_finite_positive_and_shifted() {
        let phi = 0.5_f64;
        let es = 0.00669438;
        let m = pj_msfn(phi.sin(), phi.cos(), es);
        assert!(m.is_finite());
        assert!(m > 0.0);
        // With nonzero eccentricity the result differs from plain cos(phi).
        assert!((m - phi.cos()).abs() > 0.0);
    }
}