oxiproj-core 0.1.2

//! 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 {
    /// Semi-major axis, in meters.
    pub a: f64,
    /// Semi-minor axis, in meters.
    pub b: f64,
    /// Reciprocal of the semi-major axis (`1/a`).
    pub ra: f64,
    /// Reciprocal of the semi-minor axis (`1/b`).
    pub rb: f64,
    /// Angular eccentricity, `asin(e)`, in radians.
    pub alpha: f64,
    /// First eccentricity, `sqrt(es)`.
    pub e: f64,
    /// First eccentricity squared.
    pub es: f64,
    /// Second eccentricity, `tan(alpha)`.
    pub e2: f64,
    /// Second eccentricity squared, `e2 * e2`.
    pub e2s: f64,
    /// Third eccentricity.
    pub e3: f64,
    /// Third eccentricity squared, `e3 * e3`.
    pub e3s: f64,
    /// `1 - es`.
    pub one_es: f64,
    /// Reciprocal of [`Self::one_es`].
    pub rone_es: f64,
    /// Flattening, `(a - b) / a`.
    pub f: f64,
    /// Second flattening.
    pub f2: f64,
    /// Third flattening.
    pub n: f64,
    /// Inverse flattening, `1 / f` (trig-derived; see
    /// [`Ellipsoid::precise_inv_flattening`] for a fidelity-preserving
    /// accessor suitable for WKT/PROJJSON round-trips).
    pub rf: f64,
    /// Reciprocal of [`Self::f2`].
    pub rf2: f64,
    /// Reciprocal of [`Self::n`].
    pub rn: f64,
    /// The exact inverse-flattening value the ellipsoid was constructed
    /// from via [`from_a_rf`] (or [`named`] for table entries defined by
    /// `rf=`), if any. `None` when the ellipsoid was built from `a`/`b`,
    /// `a`/`f`, or `a`/`es` directly. See
    /// [`Ellipsoid::precise_inv_flattening`].
    literal_rf: Option<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,
            literal_rf: None,
        }
    }
}

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)
    }

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

    /// Returns the most precise available inverse flattening (`1/f`).
    ///
    /// When this ellipsoid was constructed via [`from_a_rf`] (directly, or
    /// indirectly through [`named`] for table entries defined by an
    /// `rf=<value>` shape parameter, or through the proj-string/WKT parsers
    /// which also bottom out in `from_a_rf`), the exact caller-supplied
    /// literal is stored and returned unchanged here. This lets WKT2/WKT1/
    /// PROJJSON exporters print e.g. `298.257223563` for WGS84 byte-for-byte
    /// instead of a value corrupted by the `asin`/`cos` round-trip that
    /// [`calc_ellipsoid_params`] performs internally to populate
    /// [`Ellipsoid::rf`] (which remains unchanged, for numerical-code
    /// backward compatibility).
    ///
    /// Otherwise (ellipsoid built from `a`/`b`, `a`/`f`, or `a`/`es`
    /// directly, so there is no original literal inverse-flattening to
    /// preserve), the value is derived straight from `es` via the identity
    /// `f = 1 - sqrt(1 - es)`, `rf = 1/f`. This avoids the extra
    /// trigonometric round-trip (`sqrt` -> `asin` -> `cos`) used internally
    /// by [`calc_ellipsoid_params`] and is therefore numerically tighter
    /// than the [`Ellipsoid::rf`] field, though it cannot in general recover
    /// a literal value bit-for-bit once it has been folded into `es`.
    ///
    /// A perfect sphere (`f == 0`) returns [`f64::INFINITY`], matching
    /// [`Ellipsoid::rf`]'s convention.
    pub fn precise_inv_flattening(&self) -> f64 {
        if let Some(rf) = self.literal_rf {
            return rf;
        }
        if self.f == 0.0 {
            return f64::INFINITY;
        }
        let f = 1.0 - (1.0 - self.es).sqrt();
        if f != 0.0 {
            1.0 / f
        } else {
            f64::INFINITY
        }
    }
}

/// Ported from src/ellps.cpp (`pj_calc_ellipsoid_params`, together with the
/// validation performed by `ellps_size`/`ellps_shape` in `src/ell_set.cpp`).
///
/// Given the semi-major axis `a` and eccentricity-squared `es`, computes all
/// derived ellipsoid parameters. This is the single choke point every public
/// constructor in this module (`from_a_es`, `from_a_b`, `from_a_rf`,
/// `from_a_f`, `sphere`, `named`, and [`try_new`]) funnels through, so the
/// validation performed here applies uniformly to all of them.
///
/// # Errors
///
/// Returns [`ProjError::IllegalArgValue`] (PROJ error code 1027, the same
/// class real PROJ raises, e.g. `pj_init_ctx: Must specify ellipsoid or
/// sphere` for `+a=0`) when:
/// - `a` is not finite or `a <= 0.0` (PROJ: `ellps_size`, "Invalid value for
///   major axis"), or
/// - `es` is not finite, `es < 0.0`, or `es >= 1.0` (PROJ: `ellps_shape`,
///   "Invalid value for es. Should be in \[0,1\[ range").
pub fn calc_ellipsoid_params(a: f64, es: f64) -> ProjResult<Ellipsoid> {
    if !(a.is_finite() && a > 0.0) {
        return Err(ProjError::IllegalArgValue);
    }
    if !(es.is_finite() && (0.0..1.0).contains(&es)) {
        return Err(ProjError::IllegalArgValue);
    }

    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.
///
/// See [`calc_ellipsoid_params`] for the validation rules applied to `a` and
/// `es`.
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.
///
/// # Errors
///
/// Returns [`ProjError::IllegalArgValue`] when `b` is not finite or not
/// strictly positive (PROJ: `ellps_shape` case `b`, "Invalid value for b.
/// Should be > 0"), or when the derived eccentricity-squared / `a` fail the
/// checks in [`calc_ellipsoid_params`] (e.g. `a <= 0`, or `b` greater than
/// `a` yielding a negative derived `es`).
pub fn from_a_b(a: f64, b: f64) -> ProjResult<Ellipsoid> {
    if !(b.is_finite() && b > 0.0) {
        return Err(ProjError::IllegalArgValue);
    }
    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.
///
/// The exact `rf` supplied here is preserved and can be recovered losslessly
/// via [`Ellipsoid::precise_inv_flattening`], even though the [`Ellipsoid::rf`]
/// field itself is recomputed through [`calc_ellipsoid_params`]'s internal
/// trigonometric round-trip and may differ in the last few ULPs.
///
/// # Errors
///
/// Returns [`ProjError::IllegalArgValue`] when `rf` is not finite or
/// `rf <= 0.0` (PROJ: `ellps_shape` case `rf`, "Invalid value for rf. Should
/// be > 0"), or when `a`/the derived `es` fail the checks in
/// [`calc_ellipsoid_params`].
pub fn from_a_rf(a: f64, rf: f64) -> ProjResult<Ellipsoid> {
    if !(rf.is_finite() && rf > 0.0) {
        return Err(ProjError::IllegalArgValue);
    }
    let f = 1.0 / rf;
    let es = 2.0 * f - f * f;
    let mut el = calc_ellipsoid_params(a, es)?;
    el.literal_rf = Some(rf);
    Ok(el)
}

/// Ported from src/ellps.cpp. Builds from semi-major axis and flattening.
///
/// # Errors
///
/// Returns [`ProjError::IllegalArgValue`] when `f` is not finite or
/// negative (PROJ: `ellps_shape` case `f`, "Invalid value for f. Should be
/// non-negative"), or when `a`/the derived `es` fail the checks in
/// [`calc_ellipsoid_params`].
pub fn from_a_f(a: f64, f: f64) -> ProjResult<Ellipsoid> {
    if !(f.is_finite() && 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`.
///
/// # Errors
///
/// Returns [`ProjError::IllegalArgValue`] when `r` is not finite or
/// `r <= 0.0`, per [`calc_ellipsoid_params`].
pub fn sphere(r: f64) -> ProjResult<Ellipsoid> {
    from_a_es(r, 0.0)
}

/// Validating constructor: builds an [`Ellipsoid`] from a semi-major axis
/// and eccentricity-squared, rejecting invalid parameters up front.
///
/// This is a documented alias of [`calc_ellipsoid_params`], provided as the
/// canonical validating entry point for callers (such as proj-string / CRS
/// parameter parsing in downstream crates) that want a clearly-named
/// fallible constructor rather than depending on the free-function name
/// `calc_ellipsoid_params`. See [`calc_ellipsoid_params`] for the exact
/// validation rules (`a > 0` and finite; `es` in `[0, 1)` and finite) and
/// the PROJ error-code mapping ([`ProjError::IllegalArgValue`], code 1027).
///
/// # Errors
///
/// See [`calc_ellipsoid_params`].
pub fn try_new(a: f64, es: f64) -> ProjResult<Ellipsoid> {
    calc_ellipsoid_params(a, es)
}

/// 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 from_a_rf_rejects_non_finite() {
        assert!(from_a_rf(6378137.0, f64::INFINITY).is_err());
        assert!(from_a_rf(6378137.0, f64::NAN).is_err());
    }

    #[test]
    fn from_a_f_rejects_non_finite() {
        assert!(from_a_f(6378137.0, f64::INFINITY).is_err());
        assert!(from_a_f(6378137.0, f64::NAN).is_err());
    }

    #[test]
    fn from_a_b_rejects_non_positive_or_non_finite_b() {
        assert!(from_a_b(6378137.0, 0.0).is_err());
        assert!(from_a_b(6378137.0, -6356752.0).is_err());
        assert!(from_a_b(6378137.0, f64::INFINITY).is_err());
        assert!(from_a_b(6378137.0, f64::NAN).is_err());
    }

    #[test]
    fn calc_ellipsoid_params_rejects_non_positive_a() {
        // PROJ: `+proj=longlat +a=0 +b=0` fails with error 1027
        // ("Must specify ellipsoid or sphere").
        assert_eq!(
            calc_ellipsoid_params(0.0, 0.0).unwrap_err(),
            ProjError::IllegalArgValue
        );
        assert!(calc_ellipsoid_params(-6378137.0, 0.0066943799901413165).is_err());
        assert!(calc_ellipsoid_params(f64::INFINITY, 0.0).is_err());
        assert!(calc_ellipsoid_params(f64::NAN, 0.0).is_err());
    }

    #[test]
    fn calc_ellipsoid_params_rejects_invalid_es() {
        assert!(calc_ellipsoid_params(6378137.0, 1.0).is_err());
        assert!(calc_ellipsoid_params(6378137.0, 1.5).is_err());
        assert!(calc_ellipsoid_params(6378137.0, -0.1).is_err());
        assert!(calc_ellipsoid_params(6378137.0, f64::NAN).is_err());
        assert!(calc_ellipsoid_params(6378137.0, f64::INFINITY).is_err());
    }

    #[test]
    fn from_a_b_rejects_negated_axes() {
        // Both axes negated (ratio preserved, but `a` is negative): must be
        // rejected, matching real PROJ's rejection of `+a=0 +b=0` and
        // negative axis values.
        assert!(from_a_b(-6378137.0, -6356752.314245179).is_err());
    }

    #[test]
    fn try_new_matches_calc_params() {
        let a = try_new(6378137.0, 0.0066943799901413165).expect("try_new");
        let b = calc_ellipsoid_params(6378137.0, 0.0066943799901413165).expect("calc_params");
        assert_eq!(a, b);
        assert!(try_new(0.0, 0.0).is_err());
        assert_eq!(
            Ellipsoid::try_new(0.0, 0.0).unwrap_err(),
            ProjError::IllegalArgValue
        );
    }

    #[test]
    fn precise_inv_flattening_preserves_literal_wgs84() {
        // The literal 298.257223563 must round-trip exactly, unlike the
        // trig-derived `rf` field which loses precision (298.25722356300156).
        let e = from_a_rf(6378137.0, 298.257223563).expect("from_a_rf");
        assert_eq!(e.precise_inv_flattening(), 298.257223563);
        assert_ne!(
            e.rf, 298.257223563,
            "rf field is expected to be trig-corrupted"
        );
    }

    #[test]
    fn precise_inv_flattening_preserves_literal_via_named() {
        // named("WGS84") goes through the `rf=` shape parameter, so the
        // literal must be preserved end-to-end.
        let e = Ellipsoid::named("WGS84").expect("WGS84");
        assert_eq!(e.precise_inv_flattening(), 298.257223563);

        let g = Ellipsoid::named("GRS80").expect("GRS80");
        assert_eq!(g.precise_inv_flattening(), 298.257222101);
    }

    #[test]
    fn precise_inv_flattening_preserves_literal_international1924() {
        // Hand-written WKT-style ellipsoid definition using an integer rf.
        let e = from_a_rf(6378388.0, 297.0).expect("from_a_rf");
        assert_eq!(e.precise_inv_flattening(), 297.0);
    }

    #[test]
    fn precise_inv_flattening_no_literal_falls_back_to_direct_formula() {
        // Built from a/b: no literal rf to preserve, but the derived value
        // must still be finite and close to the trig-derived `rf` field.
        let e = from_a_b(6378137.0, 6356752.314245179).expect("from_a_b");
        let precise = e.precise_inv_flattening();
        assert!(precise.is_finite());
        assert!((precise - e.rf).abs() < 1e-6);
    }

    #[test]
    fn precise_inv_flattening_sphere_is_infinite() {
        let s = sphere(6370997.0).expect("sphere");
        assert!(s.precise_inv_flattening().is_infinite());
    }

    #[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());
    }
}