oxiproj-core 0.1.0

//! Ellipsoid parameter computation ported from PROJ `src/ellps.cpp`.

use crate::error::{ProjError, ProjResult};

/// A fully-computed ellipsoid (all derived parameters). Ported from src/ellps.cpp
/// (`PJ` ellipsoid fields populated by `pj_calc_ellipsoid_params`).
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct Ellipsoid {
    pub a: f64,
    pub b: f64,
    pub ra: f64,
    pub rb: f64,
    pub alpha: f64,
    pub e: f64,
    pub es: f64,
    pub e2: f64,
    pub e2s: f64,
    pub e3: f64,
    pub e3s: f64,
    pub one_es: f64,
    pub rone_es: f64,
    pub f: f64,
    pub f2: f64,
    pub n: f64,
    pub rf: f64,
    pub rf2: f64,
    pub rn: f64,
}

impl Ellipsoid {
    /// All-zero ellipsoid (internal scratch before parameter calculation).
    fn zeroed() -> Ellipsoid {
        Ellipsoid {
            a: 0.0,
            b: 0.0,
            ra: 0.0,
            rb: 0.0,
            alpha: 0.0,
            e: 0.0,
            es: 0.0,
            e2: 0.0,
            e2s: 0.0,
            e3: 0.0,
            e3s: 0.0,
            one_es: 0.0,
            rone_es: 0.0,
            f: 0.0,
            f2: 0.0,
            n: 0.0,
            rf: 0.0,
            rf2: 0.0,
            rn: 0.0,
        }
    }
}

impl Ellipsoid {
    /// Builder: see [`calc_ellipsoid_params`].
    pub fn calc_params(a: f64, es: f64) -> ProjResult<Ellipsoid> {
        crate::ellipsoid::calc_ellipsoid_params(a, es)
    }

    /// Builder: see [`from_a_es`].
    pub fn from_a_es(a: f64, es: f64) -> ProjResult<Ellipsoid> {
        crate::ellipsoid::from_a_es(a, es)
    }

    /// Builder: see [`from_a_b`].
    pub fn from_a_b(a: f64, b: f64) -> ProjResult<Ellipsoid> {
        crate::ellipsoid::from_a_b(a, b)
    }

    /// Builder: see [`from_a_rf`].
    pub fn from_a_rf(a: f64, rf: f64) -> ProjResult<Ellipsoid> {
        crate::ellipsoid::from_a_rf(a, rf)
    }

    /// Builder: see [`from_a_f`].
    pub fn from_a_f(a: f64, f: f64) -> ProjResult<Ellipsoid> {
        crate::ellipsoid::from_a_f(a, f)
    }

    /// Builder: see [`sphere`].
    pub fn sphere(r: f64) -> ProjResult<Ellipsoid> {
        crate::ellipsoid::sphere(r)
    }

    /// Builder: see [`named`].
    pub fn named(id: &str) -> ProjResult<Ellipsoid> {
        crate::ellipsoid::named(id)
    }
}

/// Ported from src/ellps.cpp (`pj_calc_ellipsoid_params`).
///
/// Given the semi-major axis `a` and eccentricity-squared `es`, computes all
/// derived ellipsoid parameters.
pub fn calc_ellipsoid_params(a: f64, es: f64) -> ProjResult<Ellipsoid> {
    let mut el = Ellipsoid::zeroed();
    el.a = a;
    el.es = es;

    // angular eccentricity & friends
    el.e = es.sqrt();
    el.alpha = el.e.asin(); // angular eccentricity
    el.e2 = el.alpha.tan();
    el.e2s = el.e2 * el.e2; // second eccentricity (squared)
    el.e3 = if el.alpha != 0.0 {
        el.alpha.sin() / (2.0 - el.alpha.sin() * el.alpha.sin()).sqrt()
    } else {
        0.0
    };
    el.e3s = el.e3 * el.e3;

    el.f = 1.0 - el.alpha.cos(); // flattening
    if !(el.f >= 0.0 && el.f < 1.0) {
        return Err(ProjError::IllegalArgValue);
    }
    el.rf = if el.f != 0.0 {
        1.0 / el.f
    } else {
        f64::INFINITY
    };

    el.f2 = if el.alpha.cos() != 0.0 {
        1.0 / el.alpha.cos() - 1.0
    } else {
        0.0
    };
    el.rf2 = if el.f2 != 0.0 {
        1.0 / el.f2
    } else {
        f64::INFINITY
    };

    el.n = (el.alpha / 2.0).tan().powi(2); // third flattening
    el.rn = if el.n != 0.0 {
        1.0 / el.n
    } else {
        f64::INFINITY
    };

    el.b = (1.0 - el.f) * el.a;
    el.rb = 1.0 / el.b;
    el.ra = 1.0 / el.a;

    el.one_es = 1.0 - el.es;
    if el.one_es == 0.0 {
        return Err(ProjError::IllegalArgValue);
    }
    el.rone_es = 1.0 / el.one_es;

    Ok(el)
}

/// Ported from src/ellps.cpp. Builds from semi-major axis and eccentricity-squared.
pub fn from_a_es(a: f64, es: f64) -> ProjResult<Ellipsoid> {
    calc_ellipsoid_params(a, es)
}

/// Ported from src/ellps.cpp. Builds from semi-major and semi-minor axes.
pub fn from_a_b(a: f64, b: f64) -> ProjResult<Ellipsoid> {
    let es = if b == a {
        0.0
    } else {
        let f = (a - b) / a;
        2.0 * f - f * f
    };
    calc_ellipsoid_params(a, es)
}

/// Ported from src/ellps.cpp. Builds from semi-major axis and inverse flattening.
pub fn from_a_rf(a: f64, rf: f64) -> ProjResult<Ellipsoid> {
    if rf <= 0.0 {
        return Err(ProjError::IllegalArgValue);
    }
    let f = 1.0 / rf;
    let es = 2.0 * f - f * f;
    calc_ellipsoid_params(a, es)
}

/// Ported from src/ellps.cpp. Builds from semi-major axis and flattening.
pub fn from_a_f(a: f64, f: f64) -> ProjResult<Ellipsoid> {
    if f < 0.0 {
        return Err(ProjError::IllegalArgValue);
    }
    let es = 2.0 * f - f * f;
    calc_ellipsoid_params(a, es)
}

/// Ported from src/ellps.cpp. Builds a perfect sphere of radius `r`.
pub fn sphere(r: f64) -> ProjResult<Ellipsoid> {
    from_a_es(r, 0.0)
}

/// Parses a numeric token like `"6376428."` (trailing dot allowed) into f64.
fn parse_axis_value(token: &str) -> Option<f64> {
    token.parse::<f64>().ok()
}

/// Ported from src/ellps.cpp. Builds a named ellipsoid from the PROJ table.
pub fn named(id: &str) -> ProjResult<Ellipsoid> {
    let def = match crate::ellps_table::find_ellps(id) {
        Some(d) => d,
        None => return Err(ProjError::NotEllipsoid),
    };
    // major is "a=<num>"
    let a_str = match def.major.strip_prefix("a=") {
        Some(s) => s,
        None => return Err(ProjError::IllegalArgValue),
    };
    let a = match parse_axis_value(a_str) {
        Some(v) => v,
        None => return Err(ProjError::IllegalArgValue),
    };
    // shape is "rf=<num>" or "b=<num>"
    if let Some(rf_str) = def.shape.strip_prefix("rf=") {
        let rf = match parse_axis_value(rf_str) {
            Some(v) => v,
            None => return Err(ProjError::IllegalArgValue),
        };
        from_a_rf(a, rf)
    } else if let Some(b_str) = def.shape.strip_prefix("b=") {
        let b = match parse_axis_value(b_str) {
            Some(v) => v,
            None => return Err(ProjError::IllegalArgValue),
        };
        from_a_b(a, b)
    } else {
        Err(ProjError::IllegalArgValue)
    }
}

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

    #[test]
    fn wgs84_known_values() {
        let e = named("WGS84").expect("WGS84");
        assert_eq!(e.a, 6378137.0);
        assert!((e.rf - 298.257223563).abs() < 1e-9);
        assert!((e.es - 0.0066943799901413165).abs() < 1e-15);
        assert!((e.b - 6356752.314245179).abs() < 1e-6);
        assert!((e.f - 0.0033528106647474805).abs() < 1e-15);
        assert!((e.n - 0.0016792203863836515).abs() < 1e-15);
    }

    #[test]
    fn grs80_es() {
        assert!((named("GRS80").expect("GRS80").es - 0.0066943800229007896).abs() < 1e-15);
    }

    #[test]
    fn sphere_is_degenerate() {
        let s = sphere(6370997.0).expect("sphere");
        assert_eq!(s.es, 0.0);
        assert_eq!(s.e, 0.0);
        assert_eq!(s.b, s.a);
        assert_eq!(s.f, 0.0);
        assert_eq!(s.n, 0.0);
        assert!(s.rf.is_infinite());
        assert!(s.rn.is_infinite());
    }

    #[test]
    fn clrk66_from_b() {
        let c = named("clrk66").expect("clrk66");
        assert_eq!(c.a, 6378206.4);
        assert!((c.es - 0.006768657997291).abs() < 1e-12);
    }

    #[test]
    fn unknown_ellipsoid_errors() {
        assert!(named("not_an_ellipsoid").is_err());
        assert_eq!(
            named("not_an_ellipsoid").unwrap_err(),
            ProjError::NotEllipsoid
        );
    }

    #[test]
    fn from_a_rf_matches_named_wgs84() {
        let a = from_a_rf(6378137.0, 298.257223563).expect("from_a_rf");
        let b = named("WGS84").expect("WGS84");
        assert!((a.es - b.es).abs() < 1e-15);
    }

    #[test]
    fn from_a_f_rejects_negative() {
        assert!(from_a_f(6378137.0, -0.1).is_err());
    }

    #[test]
    fn from_a_rf_rejects_non_positive() {
        assert!(from_a_rf(6378137.0, 0.0).is_err());
    }

    #[test]
    fn assoc_named_wgs84() {
        let e = Ellipsoid::named("WGS84").expect("WGS84");
        assert_eq!(e.a, 6378137.0);
        assert!((e.rf - 298.257223563).abs() < 1e-9);
    }

    #[test]
    fn assoc_sphere() {
        let s = Ellipsoid::sphere(6370997.0).expect("sphere");
        assert_eq!(s.es, 0.0);
        assert_eq!(s.b, 6370997.0);
    }

    #[test]
    fn assoc_from_a_rf_matches_named() {
        let a = Ellipsoid::from_a_rf(6378137.0, 298.257223563).expect("from_a_rf");
        let b = Ellipsoid::named("WGS84").expect("WGS84");
        assert!((a.a - b.a).abs() < 1e-9);
        assert!((a.es - b.es).abs() < 1e-9);
        assert!((a.b - b.b).abs() < 1e-9);
    }

    #[test]
    fn assoc_from_a_es() {
        let e = Ellipsoid::from_a_es(6378137.0, 0.0066943799901413165).expect("from_a_es");
        assert!(e.es.is_finite());
        assert_eq!(e.a, 6378137.0);
    }

    #[test]
    fn assoc_from_a_b() {
        let e = Ellipsoid::from_a_b(6378137.0, 6356752.314245179).expect("from_a_b");
        assert_eq!(e.a, 6378137.0);
        assert!((e.b - 6356752.314245179).abs() < 1e-6);
    }

    #[test]
    fn assoc_from_a_f() {
        let e = Ellipsoid::from_a_f(6378137.0, 0.0033528106647474805).expect("from_a_f");
        assert_eq!(e.a, 6378137.0);
        assert!(e.f >= 0.0);
    }

    #[test]
    fn assoc_calc_params() {
        let e = Ellipsoid::calc_params(6378137.0, 0.0066943799901413165).expect("calc_params");
        assert_eq!(e.a, 6378137.0);
        assert!(e.es.is_finite());
    }
}