oxiproj-engine 0.1.2

//! Single-operation builder: turns parsed parameters into a [`crate::pj::Pj`].
//!
//! Mirrors PROJ 9.8.0 `src/pj_init.cpp` generic parameter setup combined with
//! the per-projection dispatch in `oxiproj-projections`.

use crate::context::Context;
use crate::params::{ParamList, ParamView};
use crate::pj::Pj;
use oxiproj_core::{Coord, Ellipsoid, Operation, ProjError, ProjResult};

/// Evaluate a `+to_meter=` / `+vto_meter=` value, honoring PROJ's ratio syntax.
///
/// Ported from PROJ 9.8.0 `src/init.cpp`: the value is a leading number,
/// optionally followed by `/<denominator>`, in which case the two are divided
/// (e.g. `1/0.3048` -> ~3.2808). A zero denominator or a non-positive result is
/// rejected with [`ProjError::IllegalArgValue`], mirroring PROJ.
fn parse_ratio_factor(s: &str) -> ProjResult<f64> {
    let s = s.trim();
    let (num, consumed) = oxiproj_core::parse_leading_f64(s).ok_or(ProjError::IllegalArgValue)?;
    let rest = s.get(consumed..).unwrap_or("");
    let value = match rest.strip_prefix('/') {
        Some(after) => {
            let (denom, _) =
                oxiproj_core::parse_leading_f64(after).ok_or(ProjError::IllegalArgValue)?;
            if denom == 0.0 {
                return Err(ProjError::IllegalArgValue);
            }
            num / denom
        }
        None => num,
    };
    if value <= 0.0 {
        return Err(ProjError::IllegalArgValue);
    }
    Ok(value)
}

/// Recover an `axisswap` operation's canonical signed-permutation array by
/// probing its forward map, for use by the pipeline optimizer.
///
/// The `axisswap` conversion (`oxiproj-transformations`) applies a pure signed
/// axis permutation: `out[i] = sign[i] * in[axis[i]]`. Feeding it a probe with
/// four distinct, positive, exactly-representable magnitudes (`1, 2, 4, 8`)
/// lets us read back both the source axis (from `|out[i]|`) and the sign (from
/// `out[i]`'s sign) for every output slot, without downcasting to the concrete
/// type or re-parsing `+order=`/`+axis=` (which would risk drifting from the
/// operation's real behaviour).
///
/// Returns the array in PROJ `+order=` form (`sign * (source_axis + 1)`), or
/// `None` if the forward map is not a clean signed permutation of the probe
/// (which should never happen for a well-formed `axisswap`, but keeps the
/// optimizer conservative — an unrecognised mapping simply never cancels).
fn probe_axisswap_order(op: &dyn Operation) -> Option<[i8; 4]> {
    // Powers of two: distinct magnitudes, each exact in f64.
    const PROBE: [f64; 4] = [1.0, 2.0, 4.0, 8.0];
    let out = op
        .forward_4d(Coord::new(PROBE[0], PROBE[1], PROBE[2], PROBE[3]))
        .ok()?;
    let v = out.v();
    let mut order = [0i8; 4];
    let mut seen = [false; 4];
    for (i, slot) in order.iter_mut().enumerate() {
        let val = v[i];
        if !val.is_finite() || val == 0.0 {
            return None;
        }
        let mag = val.abs();
        // Match |val| against the probe magnitudes to recover the source axis.
        let axis = PROBE.iter().position(|&p| p == mag)?;
        if seen[axis] {
            // Two output slots claim the same input axis: not a permutation.
            return None;
        }
        seen[axis] = true;
        let sign: i8 = if val < 0.0 { -1 } else { 1 };
        *slot = sign * (axis as i8 + 1);
    }
    Some(order)
}

/// Decide whether a projection's AD-capable `ProjectGeneric` forward yields
/// PROJ-consistent exact distortion factors for this ellipsoid, so the engine
/// can safely activate the exact-AD Tissot path (`Pj::factors_exact`) instead
/// of the numeric finite-difference fallback.
///
/// This is a *principled, per-projection* gate — not a hardcoded name list. It
/// compares, at an interior sample point, the exact factors
/// ([`oxiproj_core::FactorsExact::from_jacobian_es`], driven by the generic
/// forward's AD Jacobian, `a = 1`, ellipsoid-corrected by `es`) against the
/// numeric factors ([`oxiproj_core::Factors::compute`], the PROJ 9.8.0
/// `src/factors.cpp` port driven by the same normalised forward `op`). It
/// returns `true` only when the two meridional/parallel scales agree to a
/// relative `1e-6`, which happens exactly when the generic forward reproduces
/// the numeric forward's geometry:
///
/// * A genuinely ellipsoidal generic forward (e.g. `merc`, `lcc`, `stere`,
///   `laea` on an ellipsoid) matches → exact path enabled.
/// * A spherical-only generic forward built on an ellipsoid yields a different
///   Jacobian than the numeric ellipsoidal forward → mismatch → numeric-gated.
/// * A projection whose `project_fwd_generic` is not AD-capable (returns
///   `UnsupportedOperation`, e.g. the exact Poder/Engsager `tmerc`) → no exact
///   sample → numeric-gated.
/// * On a sphere (`es == 0`) both paths use the identical spherical formula, so
///   every sphere projection that was enabled before stays enabled.
///
/// Samples several interior latitudes so a projection whose domain excludes one
/// candidate is still evaluated where it is defined; declines conservatively if
/// no in-domain, AD-capable sample is found (never worse than numeric).
fn ad_forward_is_proj_consistent(
    ad: &dyn oxiproj_core::ProjectGenericBox,
    op: &dyn Operation,
    es: f64,
    one_es: f64,
) -> bool {
    use oxiproj_core::autodiff::Dual1;

    // Relative tolerance: a true match agrees to ~1e-9 (finite-difference
    // truncation), while a spherical-vs-ellipsoidal mismatch differs by ~es/2
    // (~3e-3 for WGS84) — cleanly separated by 1e-6.
    const TOL: f64 = 1e-6;
    // Longitude offset from the central meridian (radians). Non-zero so the
    // parallel scale is exercised at an off-meridian point.
    const LAM_REL: f64 = 0.08;
    // Interior sample latitudes (radians): non-equatorial (so the ellipsoid
    // correction is observable — at the equator spherical and ellipsoidal
    // scales coincide), both hemispheres, away from the poles.
    const CANDIDATE_PHI: [f64; 12] = [
        0.5, -0.5, 0.7, -0.7, 0.9, -0.9, 0.3, -0.3, 1.1, -1.1, 0.2, -0.2,
    ];

    for &phi in &CANDIDATE_PHI {
        // Numeric factors from the raw normalised-space forward, mirroring
        // `Pj::factors` exactly (same closure, same es/one_es).
        let numeric =
            match oxiproj_core::Factors::compute(phi, LAM_REL, es, one_es, |lam_f, phi_f| {
                let r = op.forward_4d(Coord::new(lam_f, phi_f, 0.0, 0.0))?;
                let rv = r.v();
                Ok((rv[0], rv[1]))
            }) {
                Ok(f) => f,
                // Sample outside this projection's domain: try another latitude.
                Err(_) => continue,
            };

        // Exact AD factors at the same sample, mirroring `Pj::factors_exact`
        // (a = 1 normalised space, ellipsoid-corrected by es).
        let lam_d = Dual1::<2>::variable(LAM_REL, 0);
        let phi_d = Dual1::<2>::variable(phi, 1);
        let (x, y) = match ad.project_fwd_dual2(lam_d, phi_d) {
            Ok(v) => v,
            // Generic forward not AD-capable here: try another latitude.
            Err(_) => continue,
        };
        let exact = oxiproj_core::FactorsExact::from_jacobian_es(
            x.d[0], x.d[1], y.d[0], y.d[1], phi, 1.0, es,
        );

        let dh = (exact.meridian_scale - numeric.h).abs();
        let dk = (exact.parallel_scale - numeric.k).abs();
        return dh <= TOL * (1.0 + numeric.h.abs()) && dk <= TOL * (1.0 + numeric.k.abs());
    }

    // No in-domain, AD-capable sample: decline the exact path conservatively.
    false
}

/// Whether a projection's forward mapping ignores the ellipsoid's eccentricity,
/// i.e. is *sphere-only* — the exact condition under which PROJ forces
/// `P->es = 0` in the projection setup (`src/projections/*.cpp`) so that its
/// map-scale factors / Tissot metric are computed on the sphere.
///
/// This is a *principled, per-projection* signal — not a hardcoded name list —
/// mirroring [`ad_forward_is_proj_consistent`]: it builds the same projection a
/// second time on a sphere of identical radius (`es = 0`) and compares the two
/// forward maps at several interior, non-equatorial sample points (where any
/// ellipsoidal correction is observable — at the equator the spherical and
/// ellipsoidal forms coincide). Both operations run in the engine's normalised
/// (`a = 1`) space, so identical output means the projection's math does not
/// depend on `es`:
///
/// * A sphere-only forward (e.g. `moll`, `sinu`, `hammer`, `eck4`) built on an
///   ellipsoid produces output bit-for-bit identical to the sphere build →
///   `true` → factors use `es = 0`, matching PROJ.
/// * A genuinely ellipsoidal forward (e.g. `merc`, `tmerc`, `lcc`, `aea`,
///   `stere`, `laea`) differs from its sphere build at any off-equator sample →
///   `false` → factors keep the ellipsoid's `es`.
///
/// Declines conservatively (`false`, keep `es`) when no in-domain sample can be
/// compared, so an ellipsoidal projection is never wrongly forced onto the
/// sphere metric.
fn projection_forward_ignores_es(op_ell: &dyn Operation, op_sph: &dyn Operation) -> bool {
    // Longitude offset from the central meridian (radians), matching the sibling
    // gate; non-zero so the parallel direction is exercised off-meridian.
    const LAM_REL: f64 = 0.08;
    // Interior, non-equatorial sample latitudes (radians), both hemispheres,
    // away from the poles. The equator is excluded on purpose: there the
    // spherical and ellipsoidal forwards coincide, so it cannot discriminate.
    const CANDIDATE_PHI: [f64; 12] = [
        0.5, -0.5, 0.7, -0.7, 0.9, -0.9, 0.3, -0.3, 1.1, -1.1, 0.2, -0.2,
    ];
    let mut compared_any = false;
    for &phi in &CANDIDATE_PHI {
        let probe = Coord::new(LAM_REL, phi, 0.0, 0.0);
        let e = match op_ell.forward_4d(probe) {
            Ok(v) => v.v(),
            // Sample outside this projection's domain: try another latitude.
            Err(_) => continue,
        };
        let s = match op_sph.forward_4d(probe) {
            Ok(v) => v.v(),
            Err(_) => continue,
        };
        if !e[0].is_finite() || !e[1].is_finite() || !s[0].is_finite() || !s[1].is_finite() {
            continue;
        }
        compared_any = true;
        // The ellipsoid correction for a genuinely ellipsoidal forward is
        // O(es) ~ 3e-3 for WGS84; a sphere-only forward agrees to machine
        // precision. A tight relative tolerance cleanly separates the two.
        let tol = 1e-9 * (1.0 + e[0].abs().max(e[1].abs()));
        if (e[0] - s[0]).abs() > tol || (e[1] - s[1]).abs() > tol {
            return false;
        }
    }
    compared_any
}

/// Parse `+lon_wrap=<center>`, PROJ's output-longitude wrap center.
///
/// Ported from PROJ 9.8.0 `src/init.cpp` (`PIN->is_long_wrap_set` /
/// `PIN->long_wrap_center`): the value is a degree-valued angle (also
/// accepting DMS syntax, same convention as `+lon_0`), converted to radians.
/// A bare `+lon_wrap` (no `=value`) defaults to a `0.0` center, matching
/// PROJ's numeric-parse fallback for a missing operand. PROJ rejects centers
/// whose magnitude is not `< 10 * M_TWOPI` ("an excessive value [would]
/// perform badly when correcting longitudes around it").
///
/// Returns `Ok(None)` when `+lon_wrap` is absent (the ordinary `+over`-gated
/// wrap applies instead; see [`crate::pj::Pj`]'s `lon_wrap_center` field).
fn parse_lon_wrap(params: &ParamList) -> ProjResult<Option<f64>> {
    let raw = match params.get_str("lon_wrap") {
        Some(s) => s,
        None => return Ok(None),
    };
    let center = if raw.trim().is_empty() {
        0.0
    } else {
        oxiproj_core::dmstor(raw).map_err(|_| ProjError::IllegalArgValue)?
    };
    // Written as a positive `is_finite` guard plus `>=` (rather than PROJ's
    // literal `!(fabs(center) < 10 * M_TWOPI)`) to avoid comparing on a
    // partially-ordered type with a negated operator, while preserving the
    // same NaN-rejecting semantics PROJ's comment calls out explicitly
    // ("written this way to error on long_wrap_center '=' NaN").
    if !center.is_finite() || center.abs() >= 10.0 * oxiproj_core::M_TWOPI {
        return Err(ProjError::IllegalArgValue);
    }
    Ok(Some(center))
}

/// Build a single (non-pipeline) operation from its name and parameters.
///
/// Parses the generic PROJ parameters (`lon_0`, `lat_0`, `x_0`, `k_0`, units,
/// `pm`, ...), then delegates the projection-specific math to
/// [`oxiproj_projections::build`] and applies any overrides it returns.
pub fn build_single_op(
    name: &str,
    params: &ParamList,
    ellipsoid: Ellipsoid,
    context: &Context,
) -> ProjResult<Pj> {
    // --- Generic parameters ---
    let lam0 = params.get_dms("lon_0").unwrap_or(0.0);
    let phi0 = params.get_dms("lat_0").unwrap_or(0.0);
    let x0 = params.get_f64("x_0").unwrap_or(0.0);
    let y0 = params.get_f64("y_0").unwrap_or(0.0);
    let z0 = params.get_f64("z_0").unwrap_or(0.0);
    let k0 = params
        .get_f64("k_0")
        .or_else(|| params.get_f64("k"))
        .unwrap_or(1.0);

    // Horizontal unit scaling. `+to_meter` supports PROJ's ratio syntax
    // (`1/0.3048`); an explicit factor takes precedence over `+units`.
    let (to_meter, fr_meter) = match params.get_str("to_meter") {
        Some(s) if !s.trim().is_empty() => {
            let v = parse_ratio_factor(s)?;
            (v, 1.0 / v)
        }
        _ => match params.get_str("units") {
            Some(u) => match oxiproj_core::find_linear_unit(u) {
                Some(ud) => (ud.factor, 1.0 / ud.factor),
                None => (1.0, 1.0),
            },
            None => (1.0, 1.0),
        },
    };

    // Vertical unit scaling. `+vto_meter` likewise supports ratio syntax.
    let (vto_meter, vfr_meter) = match params.get_str("vto_meter") {
        Some(s) if !s.trim().is_empty() => {
            let v = parse_ratio_factor(s)?;
            (v, 1.0 / v)
        }
        _ => match params.get_str("vunits") {
            Some(u) => match oxiproj_core::find_linear_unit(u) {
                Some(ud) => (ud.factor, 1.0 / ud.factor),
                None => (1.0, 1.0),
            },
            None => (1.0, 1.0),
        },
    };

    let over = params.get_bool("over");
    let geoc = params.get_bool("geoc");
    let lon_wrap_center = parse_lon_wrap(params)?;

    let from_greenwich = match params.get_str("pm") {
        Some(pm) => oxiproj_core::prime_meridian_offset(pm)
            .or_else(|_| oxiproj_core::dmstor(pm))
            .unwrap_or(0.0),
        None => 0.0,
    };

    // --- Projection-specific build ---
    let view = ParamView(params);
    let pp = oxiproj_projections::ProjParams {
        ellipsoid: &ellipsoid,
        phi0,
        k0,
        params: &view,
    };
    match oxiproj_projections::build(name, &pp) {
        Ok(build) => {
            // Effective eccentricity for the map-scale factor / Tissot metric.
            // PROJ forces `P->es = 0` in the setup of every sphere-only
            // projection, computing its factors on the sphere even when the
            // user supplied an ellipsoid. Mirror that by detecting empirically
            // whether this projection's forward ignores `es` (see
            // [`projection_forward_ignores_es`]); if so, the factors metric uses
            // `es = 0`. Genuinely ellipsoidal projections keep `ellipsoid.es`.
            let factors_es = if ellipsoid.es != 0.0 {
                let sphere_only = Ellipsoid::sphere(ellipsoid.a)
                    .ok()
                    .and_then(|sph| {
                        let pp_sph = oxiproj_projections::ProjParams {
                            ellipsoid: &sph,
                            phi0,
                            k0,
                            params: &view,
                        };
                        oxiproj_projections::build(name, &pp_sph).ok().map(|b| {
                            projection_forward_ignores_es(&*build.operation, &*b.operation)
                        })
                    })
                    .unwrap_or(false);
                if sphere_only {
                    0.0
                } else {
                    ellipsoid.es
                }
            } else {
                0.0
            };

            // Apply overrides supplied by the projection.
            let lam0 = build.lam0_override.unwrap_or(lam0);
            let k0 = build.k0_override.unwrap_or(k0);
            let x0 = build.x0_override.unwrap_or(x0);
            let y0 = build.y0_override.unwrap_or(y0);
            let phi0 = build.phi0_override.unwrap_or(phi0);
            // Some projections force `P->over` internally (e.g. calcofi sets
            // `P->over = 1`), independent of the user proj-string.
            let over = build.over_override.unwrap_or(over);

            // Wire the AD-capable projection (registered by the constructor on
            // `ProjBuild::ad_proj`) into the built `Pj` so `Pj::factors_exact`
            // can compute exact Tissot distortion via forward-mode automatic
            // differentiation instead of the numeric finite-difference path.
            //
            // `FactorsExact::from_jacobian_es` now applies the ellipsoid
            // radius-of-curvature correction that PROJ's `factors.cpp` applies,
            // and `Pj::factors_exact` passes `a = 1` for the normalised forward,
            // so the exact path is PROJ-consistent on the ellipsoid too — but
            // *only* for projections whose `project_fwd_generic` is genuinely
            // the projection's own (ellipsoidal) forward. Some generic forwards
            // are spherical-only, or not AD-capable at all; enabling the exact
            // path for those would return non-PROJ-consistent distortion.
            //
            // Rather than a blanket regime gate or a hardcoded name list, we
            // decide per projection *empirically*: the exact AD factors are
            // activated only when they reproduce the numeric `factors()` (the
            // PROJ port) at an interior sample point. See
            // [`ad_forward_is_proj_consistent`]. Every other projection
            // transparently uses the numeric factors via the
            // `DistortionRaster`/`factors_exact` fallback.
            let ad_proj = build.ad_proj.filter(|ad| {
                ad_forward_is_proj_consistent(
                    &**ad,
                    &*build.operation,
                    ellipsoid.es,
                    ellipsoid.one_es,
                )
            });

            Ok(Pj {
                operation: build.operation,
                ellipsoid,
                factors_es,
                lam0,
                phi0,
                x0,
                y0,
                z0,
                k0,
                to_meter,
                fr_meter,
                vto_meter,
                vfr_meter,
                from_greenwich,
                over,
                geoc,
                lon_wrap_center,
                is_latlong: build.is_latlong,
                left: build.left,
                right: build.right,
                inverted: false,
                bypass_prepare_finalize: false,
                omit_fwd: false,
                omit_inv: false,
                ad_proj,
                op_name: name.to_string(),
                // A projection is never an `axisswap` conversion.
                axisswap_order: None,
            })
        }
        Err(oxiproj_core::ProjError::InvalidOp) => {
            // Not a projection: build it as a transformation/conversion.
            // Transforms are 4D operations working in cartesian/radians space;
            // they do their own coordinate-space handling in forward_4d/inverse_4d,
            // so projection-style prepare/finalize (lam0 subtraction, a-scaling,
            // false easting/northing) must NOT be applied. Mirrors PROJ driving
            // these by P->left / P->right (src/4D_api.cpp).
            let tp = oxiproj_transformations::TransParams {
                ellipsoid: &ellipsoid,
                params: &view,
                registry: Some(context as &dyn oxiproj_transformations::GridRegistry),
            };
            let tb = oxiproj_transformations::build(name, &tp)?;
            // For `axisswap`, capture the canonical signed permutation now (by
            // probing the built operation) so the pipeline optimizer can decide
            // cancellation from the ACTUAL permutation + sign, not from the
            // wrapper fields (which are all zeroed on the transformation path).
            let axisswap_order = if name == "axisswap" {
                probe_axisswap_order(&*tb.operation)
            } else {
                None
            };
            // `cart` and `geocent` are the geodetic <-> geocentric-Cartesian
            // conversions. In PROJ they are ordinary PJ operations driven by
            // `pj_fwd`/`pj_inv`, so their generic `fwd_prepare`/`fwd_finalize`
            // applies the central-meridian subtraction and, crucially, the
            // Cartesian `to_meter`/`fr_meter` output scaling (PROJ #1053). The
            // other conversions/transformations (helmert, axisswap, unitconvert,
            // push/pop, ...) do all their own coordinate-space handling and must
            // keep bypassing prepare/finalize. So route ONLY these two through
            // the generic prepare/finalize with the parsed unit factors and
            // central meridian, leaving every other transformation untouched.
            let cartesian_finalize = matches!(name, "cart" | "geocent");
            if cartesian_finalize {
                Ok(Pj {
                    operation: tb.operation,
                    ellipsoid,
                    factors_es: ellipsoid.es,
                    // `+lon_0` is honored by PROJ's fwd_prepare for cart/geocent.
                    lam0,
                    phi0: 0.0,
                    x0: 0.0,
                    y0: 0.0,
                    z0: 0.0,
                    k0: 1.0,
                    to_meter,
                    fr_meter,
                    vto_meter,
                    vfr_meter,
                    from_greenwich,
                    over,
                    // `+geoc` is a longitude/latitude concern; cart/geocent take
                    // geodetic input directly, so the geoc flag does not apply.
                    geoc: false,
                    lon_wrap_center,
                    is_latlong: false,
                    left: tb.left,
                    right: tb.right,
                    inverted: false,
                    bypass_prepare_finalize: false,
                    omit_fwd: false,
                    omit_inv: false,
                    ad_proj: None,
                    op_name: name.to_string(),
                    axisswap_order,
                })
            } else {
                Ok(Pj {
                    operation: tb.operation,
                    ellipsoid,
                    // A transformation/conversion has no lon/lat boundary and
                    // thus no map-scale factors; the metric simply carries the
                    // ellipsoid's `es`.
                    factors_es: ellipsoid.es,
                    lam0: 0.0,
                    phi0: 0.0,
                    x0: 0.0,
                    y0: 0.0,
                    z0: 0.0,
                    k0: 1.0,
                    to_meter: 1.0,
                    fr_meter: 1.0,
                    vto_meter: 1.0,
                    vfr_meter: 1.0,
                    from_greenwich: 0.0,
                    over: false,
                    geoc: false,
                    lon_wrap_center,
                    is_latlong: false,
                    left: tb.left,
                    right: tb.right,
                    inverted: false,
                    bypass_prepare_finalize: true,
                    omit_fwd: false,
                    omit_inv: false,
                    ad_proj: None,
                    op_name: name.to_string(),
                    axisswap_order,
                })
            }
        }
        Err(e) => Err(e),
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::params::parse;
    use oxiproj_core::{Coord, DEG_TO_RAD};

    #[test]
    fn merc_known_value() {
        let ell = Ellipsoid::named("WGS84").unwrap();
        let pj = build_single_op(
            "merc",
            &parse("+proj=merc +ellps=WGS84"),
            ell,
            &Context::new(),
        )
        .unwrap();
        let out = pj
            .forward(Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0))
            .unwrap();
        let o = out.v();
        assert!((o[0] - 1335833.8895192828).abs() < 1e-6, "x got {}", o[0]);
        assert!(
            (o[1] - 7_326_837.715_045_549).abs() < 1e-6,
            "y got {}",
            o[1]
        );
    }

    #[test]
    fn to_meter_ratio_parses() {
        // PROJ ratio syntax: 1/0.3048 == 3.2808..., not 1.0 (the old lossy parse).
        assert!((parse_ratio_factor("1/0.3048").unwrap() - 3.280839895013123).abs() < 1e-12);
        assert!((parse_ratio_factor("0.3048").unwrap() - 0.3048).abs() < 1e-15);
        assert!(parse_ratio_factor("1/0").is_err(), "zero denominator");
        assert!(parse_ratio_factor("0").is_err(), "non-positive factor");
    }

    #[test]
    fn to_meter_ratio_scales_output() {
        // With +to_meter=1/0.3048, fr_meter = 0.3048, so projected output is
        // the plain-metre result scaled by 0.3048 (verified against PROJ `proj`).
        let ell = Ellipsoid::named("WGS84").unwrap();
        let pj = build_single_op(
            "merc",
            &parse("+proj=merc +ellps=WGS84 +to_meter=1/0.3048"),
            ell,
            &Context::new(),
        )
        .unwrap();
        let out = pj
            .forward(Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0))
            .unwrap();
        let o = out.v();
        assert!(
            (o[0] - 0.3048 * 1335833.8895192828).abs() < 1e-3,
            "x got {}",
            o[0]
        );
        assert!(
            (o[1] - 0.3048 * 7_326_837.715_045_549).abs() < 1e-3,
            "y got {}",
            o[1]
        );
    }

    #[test]
    fn utm_central_meridian_origin() {
        let ell = Ellipsoid::named("WGS84").unwrap();
        let pj = build_single_op(
            "utm",
            &parse("+proj=utm +zone=32 +ellps=WGS84"),
            ell,
            &Context::new(),
        )
        .unwrap();
        let out = pj
            .forward(Coord::new(9.0 * DEG_TO_RAD, 0.0, 0.0, 0.0))
            .unwrap();
        let o = out.v();
        assert!((o[0] - 500000.0).abs() < 1e-6, "x got {}", o[0]);
        assert!((o[1] - 0.0).abs() < 1e-6, "y got {}", o[1]);
    }
}