oxiproj-engine 0.1.2

//! Ellipsoid resolution from proj-string parameters.
//!
//! Ported from PROJ 9.8.0 `src/ell_set.cpp` (`pj_ell_set`).

use crate::params::ParamList;
use oxiproj_core::{Ellipsoid, ProjError, ProjResult};

/// Default semi-major axis (WGS84 / GRS80 `a`) used when no size is supplied.
const DEFAULT_A: f64 = 6378137.0;

/// Resolve an [`Ellipsoid`] from the parsed proj-string parameters.
///
/// Resolution order mirrors `pj_ell_set`:
/// 1. `+R` overrules everything (a sphere).
/// 2. A baseline from `+ellps` or `+datum` supplies `a`/`es`.
/// 3. `+a` overrides the size; shape comes from the first present of
///    `+rf`, `+f`, `+es`, `+e`, `+b` (in that order).
/// 4. Pragmatic fallbacks fill any gaps (see inline notes).
///
/// Ported from `src/ell_set.cpp`.
pub fn setup_ellipsoid(params: &ParamList) -> ProjResult<Ellipsoid> {
    // 1. +R overrules everything.
    if let Some(r) = params.get_f64("R") {
        return Ellipsoid::sphere(r);
    }

    let mut a: Option<f64> = None;
    let mut es: Option<f64> = None;
    let mut have_baseline = false;
    let mut have_shape = false;
    // Explicit flattening, when the shape came from `+rf`/`+f`/`+b`. PROJ's
    // `pj_calc_ellipsoid_params` (src/ell_set.cpp) validates the flattening
    // itself against `[0, 1)` — this catches shapes that yield a legal `es`
    // yet an illegal flattening, e.g. `+f=2` (es == 0 but f == 2). `None` when
    // the shape was given as `es`/`e` (flattening is then derived from `es`,
    // which is already range-checked below).
    let mut f_shape: Option<f64> = None;

    // 2. Baseline from +ellps, else +datum (resolved through its ellipse id).
    if let Some(name) = params.get_str("ellps") {
        let base = Ellipsoid::named(name)?;
        a = Some(base.a);
        es = Some(base.es);
        have_baseline = true;
    } else if let Some(name) = params.get_str("datum") {
        if let Some(d) = oxiproj_core::find_datum(name) {
            let base = Ellipsoid::named(d.ellipse_id)?;
            a = Some(base.a);
            es = Some(base.es);
            have_baseline = true;
            // towgs84 and nadgrids datum shifts are injected as sub-pipelines in create.rs.
        }
        // If the datum is unknown, leave the baseline unset (do not error here).
    }

    // 3a. Size override.
    if let Some(av) = params.get_f64("a") {
        a = Some(av);
    }

    // 3b. Shape override; first present wins, in this fixed order
    // (rf, f, es, e, b). Each per-parameter range check mirrors PROJ 9.8.0
    // `ellps_shape` (src/ell_set.cpp) exactly, so that setup rejects the same
    // malformed shapes at create()-time that PROJ does, with
    // `IllegalArgValue` (PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE).
    if let Some(rf) = params.get_f64("rf") {
        // PROJ: "Invalid value for rf. Should be > 0" — rejects `+rf=0`,
        // negatives, and NaN (PROJ's `HUGE_VAL == rf || rf <= 0`).
        if rf.is_nan() || rf <= 0.0 {
            return Err(ProjError::IllegalArgValue);
        }
        let f = 1.0 / rf;
        es = Some(2.0 * f - f * f);
        f_shape = Some(f);
        have_shape = true;
    }
    if !have_shape {
        if let Some(f) = params.get_f64("f") {
            // PROJ: "Invalid value for f. Should be >= 0". Note `+f=2` passes
            // this gate (f >= 0) but is later rejected by the flattening range
            // check below, exactly as PROJ's `pj_calc_ellipsoid_params` does.
            if f.is_nan() || f < 0.0 {
                return Err(ProjError::IllegalArgValue);
            }
            es = Some(2.0 * f - f * f);
            f_shape = Some(f);
            have_shape = true;
        }
    }
    if !have_shape {
        if let Some(e_s) = params.get_f64("es") {
            // PROJ: "Invalid value for es. Should be in [0,1[ range".
            // `!contains` also rejects NaN (`contains(&NaN)` is false).
            if !(0.0..1.0).contains(&e_s) {
                return Err(ProjError::IllegalArgValue);
            }
            es = Some(e_s);
            have_shape = true;
        }
    }
    if !have_shape {
        if let Some(e) = params.get_f64("e") {
            // PROJ: "Invalid value for e. Should be in [0,1[ range" — rejects
            // `+e=-0.5`, `+e=1`, and NaN (`contains(&NaN)` is false).
            if !(0.0..1.0).contains(&e) {
                return Err(ProjError::IllegalArgValue);
            }
            es = Some(e * e);
            have_shape = true;
        }
    }
    if !have_shape {
        if let Some(bv) = params.get_f64("b") {
            // PROJ: "Invalid value for b. Should be > 0" (also rejects NaN).
            if bv.is_nan() || bv <= 0.0 {
                return Err(ProjError::IllegalArgValue);
            }
            match a {
                Some(av) => {
                    let f = (av - bv) / av;
                    es = Some(2.0 * f - f * f);
                    f_shape = Some(f);
                    have_shape = true;
                }
                None => return Err(ProjError::IllegalArgValue),
            }
        }
    }

    // 4. Resolve final a/es.
    //
    // Default: a single op with no size/shape/datum and no +R is GRS80. PROJ
    // 9.8.0 (`src/init.cpp`) appends `ellps=GRS80` as the default ellipsoid
    // whenever no datum/ellps/ellipsoid info is set, so an unqualified
    // `+proj=merc` runs on GRS80 (f=1/298.257222101), not WGS84.
    let (a_final, es_final) = if a.is_none() && !have_baseline && !have_shape {
        let grs80 = Ellipsoid::named("GRS80")?;
        (grs80.a, grs80.es)
    } else if a.is_some() && !have_baseline && !have_shape {
        // If +a was given with no shape and no baseline, treat it as a sphere.
        (a.unwrap_or(DEFAULT_A), 0.0)
    } else {
        (a.unwrap_or(DEFAULT_A), es.unwrap_or(0.0))
    };

    // 4b. Flattening range check, mirroring PROJ `pj_calc_ellipsoid_params`
    //     (src/ell_set.cpp): `if (!(P->f >= 0.0 && P->f < 1.0)) error`. PROJ
    //     runs this AFTER shape resolution but BEFORE spherification. The
    //     effective flattening is the explicit one when the shape came from
    //     `+rf`/`+f`/`+b`, else it is derived from `es` (`f = 1 - sqrt(1-es)`).
    //     This is what rejects `+f=2` (a legal `es == 0` but `f == 2`), which
    //     the `Ellipsoid::from_a_es` `es`-range check alone cannot catch.
    let f_effective = match f_shape {
        Some(f) => f,
        None => 1.0 - (1.0 - es_final).sqrt(),
    };
    if !(0.0..1.0).contains(&f_effective) {
        return Err(ProjError::IllegalArgValue);
    }

    // 5. Spherification flags (+R_A, +R_V, +R_a, +R_g, +R_h, +R_lat_a,
    //    +R_lat_g, +R_C). Ported from `ellps_spherification` in
    //    `src/ell_set.cpp`. A `+R` size parameter (handled by the early return
    //    above) suppresses shape and spherification, matching PROJ.
    let (a_final, es_final) = apply_spherification(params, a_final, es_final)?;

    Ellipsoid::from_a_es(a_final, es_final)
}

/// Apply an ellipsoid spherification flag, mapping the ellipsoid onto a sphere.
///
/// Ported verbatim from `ellps_spherification` in PROJ 9.8.0 `src/ell_set.cpp`.
/// The keys are tested in PROJ's fixed precedence order and the first present
/// one wins; the result is always a sphere (`es == 0`). Returns the input
/// unchanged when no spherification key is present.
fn apply_spherification(params: &ParamList, a: f64, es: f64) -> ProjResult<(f64, f64)> {
    // Series coefficients from `src/ell_set.cpp`.
    const SIXTH: f64 = 1.0 / 6.0;
    const RA4: f64 = 17.0 / 360.0;
    const RA6: f64 = 67.0 / 3024.0;
    const RV4: f64 = 5.0 / 72.0;
    const RV6: f64 = 55.0 / 1296.0;

    let b = a * (1.0 - es).sqrt();

    // `R_lat_a`/`R_lat_g`/`R_C` share a latitude-dependent radius formula.
    let lat_radius = |phi: f64, geometric: bool| -> ProjResult<f64> {
        if phi.abs() > oxiproj_core::M_HALFPI {
            return Err(ProjError::IllegalArgValue);
        }
        let s = phi.sin();
        let t = 1.0 - es * s * s;
        if t == 0.0 {
            return Err(ProjError::IllegalArgValue);
        }
        if geometric {
            Ok(a * (1.0 - es).sqrt() / t)
        } else {
            Ok(a * (1.0 - es + t) / (2.0 * t * t.sqrt()))
        }
    };

    let new_a = if params.exists("R_A") {
        // Sphere with the same surface area as the ellipsoid.
        a * (1.0 - es * (SIXTH + es * (RA4 + es * RA6)))
    } else if params.exists("R_V") {
        // Sphere with the same volume as the ellipsoid.
        a * (1.0 - es * (SIXTH + es * (RV4 + es * RV6)))
    } else if params.exists("R_a") {
        // Arithmetic mean R = (a + b) / 2.
        (a + b) / 2.0
    } else if params.exists("R_g") {
        // Geometric mean R = sqrt(a * b).
        (a * b).sqrt()
    } else if params.exists("R_h") {
        // Harmonic mean R = 2ab / (a + b).
        if a + b == 0.0 {
            return Err(ProjError::IllegalArgValue);
        }
        (2.0 * a * b) / (a + b)
    } else if params.exists("R_lat_a") {
        // Arithmetic mean of the ellipsoid radii at the given latitude.
        let phi = params
            .get_dms("R_lat_a")
            .ok_or(ProjError::IllegalArgValue)?;
        lat_radius(phi, false)?
    } else if params.exists("R_lat_g") {
        // Geometric mean of the ellipsoid radii at the given latitude.
        let phi = params
            .get_dms("R_lat_g")
            .ok_or(ProjError::IllegalArgValue)?;
        lat_radius(phi, true)?
    } else if params.exists("R_C") {
        // Radius of the conformal sphere. PROJ's `ellps_spherification`
        // (src/ell_set.cpp, case 7) evaluates this at `P->phi0`, but the
        // spherification runs inside `pj_ellipsoid()` (src/init.cpp:566), which
        // is called BEFORE `P->phi0` is parsed from `+lat_0` (src/init.cpp:651).
        // So `P->phi0` is still 0 at that point and R_C is always taken at the
        // equator, giving R = a·√(1−es) (the geometric/semi-minor radius),
        // independent of any `+lat_0`. Match that exactly (verified against
        // `+proj=merc +R_C +ellps=WGS84 +lat_0=45`, whose x = b·λ, not the
        // conformal radius at 45°).
        lat_radius(0.0, true)?
    } else {
        // No spherification requested.
        return Ok((a, es));
    };

    if new_a <= 0.0 {
        return Err(ProjError::IllegalArgValue);
    }
    Ok((new_a, 0.0))
}

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

    #[test]
    fn wgs84_by_ellps() {
        let e = setup_ellipsoid(&parse("+ellps=WGS84")).unwrap();
        assert_eq!(e.a, 6378137.0);
        assert!((e.es - 0.0066943799901413165).abs() < 1e-15);
    }

    #[test]
    fn default_ellipsoid_is_grs80_not_wgs84() {
        // PROJ 9.8.0 (`src/init.cpp`) appends `ellps=GRS80` as the default when
        // no datum/ellps/ellipsoid info is set, so an unqualified `+proj=merc`
        // runs on GRS80, not WGS84. GRS80: a=6378137, f=1/298.257222101.
        let e = setup_ellipsoid(&parse("+proj=merc")).unwrap();
        let grs80 = Ellipsoid::named("GRS80").unwrap();
        let wgs84 = Ellipsoid::named("WGS84").unwrap();
        assert_eq!(e.a, 6378137.0);
        assert!(
            (e.es - grs80.es).abs() < 1e-18,
            "default es must equal GRS80 ({}), got {}",
            grs80.es,
            e.es
        );
        // GRS80 and WGS84 differ in es at ~3e-11; ensure we did NOT pick WGS84.
        assert!(
            (e.es - wgs84.es).abs() > 1e-12,
            "default must be GRS80, not WGS84 (es must differ)"
        );
    }

    #[test]
    fn sphere_by_r() {
        let e = setup_ellipsoid(&parse("+R=6371000")).unwrap();
        assert_eq!(e.es, 0.0);
        assert_eq!(e.a, 6371000.0);
    }

    #[test]
    fn a_and_rf() {
        let e = setup_ellipsoid(&parse("+a=6378137 +rf=298.257223563")).unwrap();
        let wgs84 = Ellipsoid::named("WGS84").unwrap();
        assert!((e.es - wgs84.es).abs() < 1e-12);
    }

    #[test]
    fn datum_resolves_ellipsoid() {
        let e = setup_ellipsoid(&parse("+datum=WGS84")).unwrap();
        assert_eq!(e.a, 6378137.0);
    }

    #[test]
    fn a_alone_is_sphere() {
        let e = setup_ellipsoid(&parse("+a=6378137")).unwrap();
        assert_eq!(e.es, 0.0);
    }

    // Spherification flags: expected radii computed from the PROJ 9.8.0
    // `ellps_spherification` formulas for WGS84 (a=6378137, es=0.0066943799901413165).
    #[test]
    fn spherify_r_a_authalic() {
        let e = setup_ellipsoid(&parse("+ellps=WGS84 +R_A")).unwrap();
        assert_eq!(e.es, 0.0, "spherified ellipsoid must be a sphere");
        assert!((e.a - 6371007.181082429).abs() < 1e-4, "R_A a={}", e.a);
        assert_eq!(e.a, e.b, "sphere: a==b");
    }

    #[test]
    fn spherify_r_v_volume() {
        let e = setup_ellipsoid(&parse("+ellps=WGS84 +R_V")).unwrap();
        assert_eq!(e.es, 0.0);
        assert!((e.a - 6371000.790396208).abs() < 1e-4, "R_V a={}", e.a);
    }

    #[test]
    fn spherify_r_a_arithmetic_mean() {
        let e = setup_ellipsoid(&parse("+ellps=WGS84 +R_a")).unwrap();
        assert_eq!(e.es, 0.0);
        assert!((e.a - 6367444.65712259).abs() < 1e-4, "R_a a={}", e.a);
    }

    #[test]
    fn spherify_r_g_geometric_mean() {
        let e = setup_ellipsoid(&parse("+ellps=WGS84 +R_g")).unwrap();
        assert_eq!(e.es, 0.0);
        assert!((e.a - 6367435.679716192).abs() < 1e-4, "R_g a={}", e.a);
    }

    #[test]
    fn spherify_r_h_harmonic_mean() {
        let e = setup_ellipsoid(&parse("+ellps=WGS84 +R_h")).unwrap();
        assert_eq!(e.es, 0.0);
        assert!((e.a - 6367426.702322451).abs() < 1e-4, "R_h a={}", e.a);
    }

    #[test]
    fn spherify_r_lat_a_at_45() {
        let e = setup_ellipsoid(&parse("+ellps=WGS84 +R_lat_a=45")).unwrap();
        assert_eq!(e.es, 0.0);
        assert!((e.a - 6378110.052870349).abs() < 1e-4, "R_lat_a a={}", e.a);
    }

    #[test]
    fn spherify_r_lat_g_at_45() {
        let e = setup_ellipsoid(&parse("+ellps=WGS84 +R_lat_g=45")).unwrap();
        assert_eq!(e.es, 0.0);
        assert!((e.a - 6378101.030201018).abs() < 1e-4, "R_lat_g a={}", e.a);
    }

    #[test]
    fn spherify_precedence_r_a_over_r_v() {
        // When multiple flags are present PROJ picks the first in its key order
        // (R_A before R_V), so the authalic radius wins.
        let e = setup_ellipsoid(&parse("+ellps=WGS84 +R_V +R_A")).unwrap();
        assert!(
            (e.a - 6371007.181082429).abs() < 1e-4,
            "precedence a={}",
            e.a
        );
    }

    #[test]
    fn spherify_ignored_when_r_given() {
        // +R fixes the radius; spherification (and shape) are ignored.
        let e = setup_ellipsoid(&parse("+R=6371000 +R_A")).unwrap();
        assert_eq!(e.a, 6371000.0);
        assert_eq!(e.es, 0.0);
    }

    #[test]
    fn spherify_r_lat_a_out_of_range_rejected() {
        assert!(setup_ellipsoid(&parse("+ellps=WGS84 +R_lat_a=91")).is_err());
    }

    // audit confirmed[11]: the engine's ellipsoid resolution must reject a
    // degenerate `a` (or a shape override that drives the derived `a`/`es`
    // out of range), matching PROJ's `pj_init_ctx: Must specify ellipsoid or
    // sphere` (error 1027). The core lane's `calc_ellipsoid_params`
    // (`oxiproj-core`) now performs this check (`a.is_finite() && a > 0.0`,
    // `es.is_finite() && (0,1)`); `setup_ellipsoid`'s single choke point —
    // the final `Ellipsoid::from_a_es(a_final, es_final)` call — inherits it
    // for every resolution path (bare `+a`, `+a`+shape override, `+b`).
    // Cross-checked live: Homebrew PROJ 9.7.0
    // `projinfo "+proj=longlat +a=0 +b=0 +no_defs +type=crs"` fails with
    // exactly that error text.

    #[test]
    fn zero_semi_major_axis_rejected() {
        // `+a=0` alone: setup.rs's step-4 resolution treats a bare `+a` (no
        // shape/baseline) as a sphere of that radius, so this reaches
        // `Ellipsoid::from_a_es(0.0, 0.0)`, which the core validation rejects.
        assert!(setup_ellipsoid(&parse("+a=0")).is_err());
    }

    #[test]
    fn zero_semi_major_axis_with_b_rejected() {
        // `+a=0 +b=0`: the degenerate all-zero ellipsoid PROJ explicitly
        // rejects (`pj_init_ctx: Must specify ellipsoid or sphere`).
        assert!(setup_ellipsoid(&parse("+a=0 +b=0")).is_err());
    }

    #[test]
    fn negative_semi_major_axis_rejected() {
        // Both axes negated (so their RATIO, and hence the derived shape, is
        // still positive) must still be rejected on `a`'s sign alone,
        // matching the live PROJ cross-check in the confirmed finding.
        assert!(setup_ellipsoid(&parse("+a=-6378137 +b=-6356752")).is_err());
    }

    #[test]
    fn negative_inverse_flattening_rejected() {
        // `+rf` negative drives the derived `es` negative, outside the
        // `calc_ellipsoid_params` `[0, 1)` range.
        assert!(setup_ellipsoid(&parse("+a=6378137 +rf=-298.257223563")).is_err());
    }

    #[test]
    fn negative_flattening_rejected() {
        // `+f` negative likewise drives `es` negative.
        assert!(setup_ellipsoid(&parse("+a=6378137 +f=-0.003")).is_err());
    }

    // Per-shape-parameter range validation, mirroring PROJ 9.8.0 `ellps_shape`
    // + `pj_calc_ellipsoid_params` (src/ell_set.cpp). Cross-checked live against
    // Homebrew PROJ 9.7.0, all failing with error 1027 (illegal arg value):
    //   proj +proj=merc +R_a +a=2 +f=2
    //   proj +proj=utm +zone=32 +ellps=GRS80 +rf=0
    //   proj +proj=utm +zone=32 +ellps=GRS80 +e=-0.5
    //   proj +proj=utm +zone=32 +ellps=GRS80 +e=1

    #[test]
    fn flattening_two_rejected() {
        // `+f=2`: es == 0 (a legal eccentricity), but the flattening 2 is out of
        // PROJ's `[0, 1)` range — caught by the `pj_calc_ellipsoid_params`
        // flattening check, NOT by the `es` range check. gie ellipsoid.gie:
        // `+proj=merc +R_a +a=2 +f=2` => failure invalid_op_illegal_arg_value.
        assert!(setup_ellipsoid(&parse("+a=2 +f=2")).is_err());
        assert!(setup_ellipsoid(&parse("+ellps=GRS80 +f=2")).is_err());
    }

    #[test]
    fn zero_inverse_flattening_rejected() {
        // `+rf=0`: PROJ rejects with "Should be > 0" (a plain sphere needs `+f=0`
        // or no shape, not `+rf=0`). gie ellipsoid.gie.
        assert!(setup_ellipsoid(&parse("+ellps=GRS80 +rf=0")).is_err());
    }

    #[test]
    fn negative_eccentricity_rejected() {
        // `+e=-0.5`: e*e == 0.25 is a legal es, but PROJ rejects e < 0
        // ("Should be in [0,1[ range"). gie ellipsoid.gie.
        assert!(setup_ellipsoid(&parse("+ellps=GRS80 +e=-0.5")).is_err());
    }

    #[test]
    fn unit_eccentricity_rejected() {
        // `+e=1` and `+es=1` are at/over the open upper bound. gie ellipsoid.gie.
        assert!(setup_ellipsoid(&parse("+ellps=GRS80 +e=1")).is_err());
        assert!(setup_ellipsoid(&parse("+ellps=GRS80 +es=1")).is_err());
    }

    #[test]
    fn negative_semiminor_axis_rejected() {
        // `+b <= 0` rejected ("Should be > 0").
        assert!(setup_ellipsoid(&parse("+a=6378137 +b=-6356752")).is_err());
        assert!(setup_ellipsoid(&parse("+a=6378137 +b=0")).is_err());
    }

    #[test]
    fn valid_shapes_still_accepted() {
        // Guard against over-eager rejection: ordinary valid shapes must pass.
        assert!(setup_ellipsoid(&parse("+a=6378137 +rf=298.257223563")).is_ok());
        assert!(setup_ellipsoid(&parse("+a=6378137 +f=0.0033528106647474805")).is_ok());
        assert!(setup_ellipsoid(&parse("+a=6378137 +es=0.006694379990141316")).is_ok());
        assert!(setup_ellipsoid(&parse("+a=6378137 +e=0.08181919104281579")).is_ok());
        assert!(setup_ellipsoid(&parse("+a=6378137 +b=6356752.314245179")).is_ok());
        // b == a is a valid sphere.
        assert!(setup_ellipsoid(&parse("+a=6378137 +b=6378137")).is_ok());
        // f == 0 is a valid sphere.
        assert!(setup_ellipsoid(&parse("+a=6378137 +f=0")).is_ok());
    }

    #[test]
    fn zero_semi_major_axis_with_valid_shape_rejected() {
        // `+a=0` combined with an otherwise-valid shape override (so
        // `have_shape` is true and step 4 does NOT take the bare-`+a`
        // sphere shortcut) must still be rejected on `a` alone.
        assert!(setup_ellipsoid(&parse("+a=0 +rf=298.257223563")).is_err());
    }
}