oxiproj-core 0.1.2

//! Guaranteed-bounds forward map projections and datum shifts (T7.4).
//!
//! Built on the outward-rounded [`Interval`] type, these transforms take
//! coordinate *intervals* and return boxes ([`IntervalCoord`] / [`IntervalXyz`])
//! that provably enclose the exact transformed coordinate. Running them on
//! degenerate point-intervals yields a certified box around a single point.
//!
//! Every transform is paired with an ordinary `f64` forward (`*_f64`) computing
//! the identical formula; by the enclosure property of [`Interval`], the box
//! from point-inputs always contains the `f64` result. The `f64` paths are
//! validated against the reference PROJ 9.8.0 implementation.
//!
//! # Coverage
//!
//! - [`MercParams`] — Mercator (ellipsoidal), closed form; matches PROJ
//!   `merc_e_forward` to rounding.
//! - [`TmercParams`] — Transverse Mercator, the Evenden/Snyder *approximate*
//!   Gauss–Krüger series (PROJ `+approx`). The meridional arc uses the classic
//!   6th-order-in-`e` series, so agreement with PROJ's `n`-based expansion is at
//!   the sub-millimetre level for `|f| <= 1/150` — see [`TmercParams::forward`].
//! - [`Helmert7`] — 7-parameter Helmert datum shift, both the approximate
//!   (linearised) and exact rotation matrices, in either sign convention;
//!   matches PROJ `helmert`.

use crate::interval::{Interval, IntervalCoord};
use crate::math_compat as mc;

/// Arc-seconds to radians (`π / 180 / 3600`), matching PROJ `ARCSEC_TO_RAD`.
pub const ARCSEC_TO_RAD: f64 = core::f64::consts::PI / 180.0 / 3600.0;

#[inline]
fn pt(x: f64) -> Interval {
    Interval::point(x)
}

// Gauss–Krüger series coefficients, identical to PROJ tmerc.cpp `FC1..FC8`.
// PROJ writes them as repeating decimals; the exact rational forms below round
// to the identical `f64` bit patterns while avoiding excessive-precision noise.
const FC1: f64 = 1.0;
const FC2: f64 = 0.5;
const FC3: f64 = 1.0 / 6.0;
const FC4: f64 = 1.0 / 12.0;
const FC5: f64 = 0.05;
const FC6: f64 = 1.0 / 30.0;
const FC7: f64 = 1.0 / 42.0;
const FC8: f64 = 1.0 / 56.0;

// ---- f64 helpers shared by the reference paths ----

#[inline]
fn asinh_f64(t: f64) -> f64 {
    mc::ln(t + mc::sqrt(t * t + 1.0))
}

#[inline]
fn atanh_f64(u: f64) -> f64 {
    0.5 * mc::ln((1.0 + u) / (1.0 - u))
}

/// Classic Snyder meridional-arc series `M/a` (dimensionless), 6th order in the
/// first eccentricity `e` (Snyder, *Map Projections — A Working Manual*, eq.
/// 3-21). Used by the transverse-Mercator forward.
#[inline]
fn meridional_arc_f64(phi: f64, es: f64) -> f64 {
    let e2 = es;
    let e4 = e2 * e2;
    let e6 = e4 * e2;
    let c0 = 1.0 - e2 / 4.0 - 3.0 * e4 / 64.0 - 5.0 * e6 / 256.0;
    let c2 = 3.0 * e2 / 8.0 + 3.0 * e4 / 32.0 + 45.0 * e6 / 1024.0;
    let c4 = 15.0 * e4 / 256.0 + 45.0 * e6 / 1024.0;
    let c6 = 35.0 * e6 / 3072.0;
    c0 * phi - c2 * mc::sin(2.0 * phi) + c4 * mc::sin(4.0 * phi) - c6 * mc::sin(6.0 * phi)
}

/// Interval meridional-arc series `M/a`; the interval mirror of
/// [`meridional_arc_f64`].
#[inline]
fn meridional_arc_iv(phi: Interval, es: f64) -> Interval {
    let e2 = pt(es);
    let e4 = e2 * e2;
    let e6 = e4 * e2;
    let c0 = pt(1.0) - e2 * pt(0.25) - e4 * pt(3.0 / 64.0) - e6 * pt(5.0 / 256.0);
    let c2 = e2 * pt(3.0 / 8.0) + e4 * pt(3.0 / 32.0) + e6 * pt(45.0 / 1024.0);
    let c4 = e4 * pt(15.0 / 256.0) + e6 * pt(45.0 / 1024.0);
    let c6 = e6 * pt(35.0 / 3072.0);
    c0 * phi - c2 * (pt(2.0) * phi).sin() + c4 * (pt(4.0) * phi).sin() - c6 * (pt(6.0) * phi).sin()
}

/// Parameters for the ellipsoidal Mercator forward projection.
///
/// Angles are in radians; `lam0` (central meridian) is subtracted from the
/// input longitude. Output is in the units of `a` (metres for a metric `a`),
/// with false easting/northing `x0`/`y0` added, exactly as PROJ.
#[derive(Clone, Copy, Debug, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct MercParams {
    /// Semi-major axis.
    pub a: f64,
    /// First eccentricity `e = sqrt(es)`.
    pub e: f64,
    /// Scale factor at the equator.
    pub k0: f64,
    /// False easting.
    pub x0: f64,
    /// False northing.
    pub y0: f64,
    /// Central meridian (radians).
    pub lam0: f64,
}

impl MercParams {
    /// Ordinary `f64` ellipsoidal Mercator forward, matching PROJ
    /// `merc_e_forward` (plus the framework's `a`-scaling and false origin).
    ///
    /// `x = a·k0·(λ − λ0) + x0`,
    /// `y = a·k0·(asinh(tanφ) − e·atanh(e·sinφ)) + y0`.
    pub fn forward_f64(&self, lam: f64, phi: f64) -> (f64, f64) {
        let dlam = lam - self.lam0;
        let sphi = mc::sin(phi);
        let cphi = mc::cos(phi);
        let iso = asinh_f64(sphi / cphi) - self.e * atanh_f64(self.e * sphi);
        let ak0 = self.a * self.k0;
        (ak0 * dlam + self.x0, ak0 * iso + self.y0)
    }

    /// Guaranteed-bounds Mercator forward.
    ///
    /// Enclosure argument: the formula is a branch-free composition of the
    /// enclosure primitives (`sin`, `cos`, `/`, `asinh`, `atanh`, `+`, `*`), so
    /// the result box provably contains the true projected point whenever the
    /// input intervals contain the true `(λ, φ)`. The latitude interval must
    /// stay strictly inside `(−π/2, π/2)` for a finite box (otherwise `cosφ`
    /// straddles 0 and the isometric latitude is unbounded — still a valid, if
    /// useless, enclosure).
    pub fn forward(&self, lam: Interval, phi: Interval) -> IntervalCoord {
        let dlam = lam - pt(self.lam0);
        let sphi = phi.sin();
        let cphi = phi.cos();
        let e = pt(self.e);
        let iso = (sphi / cphi).asinh() - e * (e * sphi).atanh();
        let ak0 = pt(self.a) * pt(self.k0);
        let x = ak0 * dlam + pt(self.x0);
        let y = ak0 * iso + pt(self.y0);
        IntervalCoord::new(x, y)
    }
}

/// Parameters for the ellipsoidal (approximate/Gauss–Krüger) Transverse
/// Mercator forward projection.
///
/// Angles are radians. `lam0`/`phi0` are the central meridian and latitude of
/// origin; output is scaled by `a` with false easting/northing added.
#[derive(Clone, Copy, Debug, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct TmercParams {
    /// Semi-major axis.
    pub a: f64,
    /// Eccentricity squared `es = e²`.
    pub es: f64,
    /// Scale factor on the central meridian.
    pub k0: f64,
    /// False easting.
    pub x0: f64,
    /// False northing.
    pub y0: f64,
    /// Central meridian (radians).
    pub lam0: f64,
    /// Latitude of origin (radians).
    pub phi0: f64,
}

impl TmercParams {
    /// Ordinary `f64` transverse-Mercator forward (Evenden/Snyder approximate
    /// series), replicating PROJ `approx_e_fwd` except that the meridional arc
    /// uses the classic 6th-order-in-`e` series (`meridional_arc_f64`).
    pub fn forward_f64(&self, lam: f64, phi: f64) -> (f64, f64) {
        let esp = self.es / (1.0 - self.es);
        let dlam = lam - self.lam0;
        let sphi = mc::sin(phi);
        let cphi = mc::cos(phi);
        let tt = if mc::abs(cphi) > 1e-10 {
            sphi / cphi
        } else {
            0.0
        };
        let t = tt * tt;
        let al0 = cphi * dlam;
        let als = al0 * al0;
        let al = al0 / mc::sqrt(1.0 - self.es * sphi * sphi);
        let n = esp * cphi * cphi;
        let x_ser = FC1
            + FC3
                * als
                * (1.0 - t
                    + n
                    + FC5
                        * als
                        * (5.0
                            + t * (t - 18.0)
                            + n * (14.0 - 58.0 * t)
                            + FC7 * als * (61.0 + t * (t * (179.0 - t) - 479.0))));
        let m = meridional_arc_f64(phi, self.es) - meridional_arc_f64(self.phi0, self.es);
        let y_ser = 1.0
            + FC4
                * als
                * (5.0 - t
                    + n * (9.0 + 4.0 * n)
                    + FC6
                        * als
                        * (61.0
                            + t * (t - 58.0)
                            + n * (270.0 - 330.0 * t)
                            + FC8 * als * (1385.0 + t * (t * (543.0 - t) - 3111.0))));
        let x_norm = self.k0 * al * x_ser;
        let y_norm = self.k0 * (m + sphi * al * dlam * FC2 * y_ser);
        (self.a * x_norm + self.x0, self.a * y_norm + self.y0)
    }

    /// Guaranteed-bounds transverse-Mercator forward.
    ///
    /// Enclosure argument: the transform is the same branch-free polynomial /
    /// meridional-series composition as [`TmercParams::forward_f64`], evaluated
    /// with enclosure primitives, so the box encloses the exact value of *that*
    /// series. Because the underlying model is the truncated Gauss–Krüger /
    /// Snyder series (not an exact closed form), the box brackets the series
    /// value rigorously; the residual series-truncation error versus the true
    /// ellipsoidal transverse Mercator is sub-millimetre for `|f| <= 1/150`
    /// (documented, matching PROJ `+approx` to that level). The latitude
    /// interval must avoid the poles (`cosφ > 0`) for a finite box.
    pub fn forward(&self, lam: Interval, phi: Interval) -> IntervalCoord {
        let esp = self.es / (1.0 - self.es);
        let dlam = lam - pt(self.lam0);
        let sphi = phi.sin();
        let cphi = phi.cos();
        let t_tan = sphi / cphi;
        let t = t_tan.powi(2);
        let al0 = cphi * dlam;
        let als = al0.powi(2);
        let denom = (pt(1.0) - pt(self.es) * sphi.powi(2)).sqrt();
        let al = al0 / denom;
        let n = pt(esp) * cphi.powi(2);
        let x_ser = pt(FC1)
            + pt(FC3)
                * als
                * ((pt(1.0) - t + n)
                    + pt(FC5)
                        * als
                        * ((pt(5.0) + t * (t - pt(18.0)) + n * (pt(14.0) - pt(58.0) * t))
                            + pt(FC7) * als * (pt(61.0) + t * (t * (pt(179.0) - t) - pt(479.0)))));
        let m = meridional_arc_iv(phi, self.es) - meridional_arc_iv(pt(self.phi0), self.es);
        let y_ser = pt(1.0)
            + pt(FC4)
                * als
                * ((pt(5.0) - t + n * (pt(9.0) + pt(4.0) * n))
                    + pt(FC6)
                        * als
                        * ((pt(61.0) + t * (t - pt(58.0)) + n * (pt(270.0) - pt(330.0) * t))
                            + pt(FC8)
                                * als
                                * (pt(1385.0) + t * (t * (pt(543.0) - t) - pt(3111.0)))));
        let x_norm = pt(self.k0) * al * x_ser;
        let y_norm = pt(self.k0) * (m + sphi * al * dlam * pt(FC2) * y_ser);
        let a = pt(self.a);
        IntervalCoord::new(a * x_norm + pt(self.x0), a * y_norm + pt(self.y0))
    }
}

/// A 3-D coordinate with guaranteed-enclosure [`Interval`] components; the
/// output of the Helmert forward.
#[derive(Clone, Copy, Debug, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct IntervalXyz {
    /// `X` (or easting) enclosure.
    pub x: Interval,
    /// `Y` (or northing) enclosure.
    pub y: Interval,
    /// `Z` (or up) enclosure.
    pub z: Interval,
}

impl IntervalXyz {
    /// Construct from three component intervals.
    #[inline]
    pub const fn new(x: Interval, y: Interval, z: Interval) -> IntervalXyz {
        IntervalXyz { x, y, z }
    }

    /// Does this box contain the point `(x, y, z)`?
    #[inline]
    pub fn contains(self, x: f64, y: f64, z: f64) -> bool {
        self.x.contains(x) && self.y.contains(y) && self.z.contains(z)
    }
}

/// Sign convention for the Helmert rotation terms (see PROJ `helmert`).
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub enum Convention {
    /// Position-vector convention (rotate the point). Used by e.g. `+towgs84`.
    PositionVector,
    /// Coordinate-frame convention (rotate the axes). Transpose of the above.
    CoordinateFrame,
}

/// Rotation-matrix construction mode (see PROJ `build_rot_matrix`).
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub enum Rotation {
    /// Small-angle linearised matrix (`cos≈1`, `sin≈angle`); PROJ default.
    Approximate,
    /// Full trigonometric rotation matrix (PROJ `+exact`).
    Exact,
}

/// 7-parameter Helmert datum-shift transformation (translation, rotation,
/// scale), matching PROJ `helmert_forward_3d` with `refp = 0`.
#[derive(Clone, Copy, Debug, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct Helmert7 {
    /// Translation `X` (metres).
    pub dx: f64,
    /// Translation `Y` (metres).
    pub dy: f64,
    /// Translation `Z` (metres).
    pub dz: f64,
    /// Rotation about `X` (radians).
    pub rx: f64,
    /// Rotation about `Y` (radians).
    pub ry: f64,
    /// Rotation about `Z` (radians).
    pub rz: f64,
    /// Scale (parts per million).
    pub ds: f64,
    /// Sign convention for the rotation terms.
    pub convention: Convention,
    /// Rotation-matrix construction mode.
    pub rotation: Rotation,
}

impl Helmert7 {
    /// Build from published parameters: translations in metres, rotations in
    /// arc-seconds, scale in ppm.
    #[allow(clippy::too_many_arguments)]
    pub fn from_arcsec_ppm(
        dx: f64,
        dy: f64,
        dz: f64,
        rx_arcsec: f64,
        ry_arcsec: f64,
        rz_arcsec: f64,
        ds_ppm: f64,
        convention: Convention,
        rotation: Rotation,
    ) -> Helmert7 {
        Helmert7 {
            dx,
            dy,
            dz,
            rx: rx_arcsec * ARCSEC_TO_RAD,
            ry: ry_arcsec * ARCSEC_TO_RAD,
            rz: rz_arcsec * ARCSEC_TO_RAD,
            ds: ds_ppm,
            convention,
            rotation,
        }
    }

    /// Rotation matrix in `f64`, exactly following PROJ `build_rot_matrix`
    /// (coordinate-frame equations, transposed for the position-vector
    /// convention).
    fn rot_f64(&self) -> [[f64; 3]; 3] {
        let (f, t, p) = (self.rx, self.ry, self.rz);
        let mut r = match self.rotation {
            Rotation::Exact => {
                let (cf, sf) = (mc::cos(f), mc::sin(f));
                let (ct, st) = (mc::cos(t), mc::sin(t));
                let (cp, sp) = (mc::cos(p), mc::sin(p));
                [
                    [ct * cp, cf * sp + sf * st * cp, sf * sp - cf * st * cp],
                    [-ct * sp, cf * cp - sf * st * sp, sf * cp + cf * st * sp],
                    [st, -sf * ct, cf * ct],
                ]
            }
            Rotation::Approximate => [[1.0, p, -t], [-p, 1.0, f], [t, -f, 1.0]],
        };
        if self.convention == Convention::PositionVector {
            transpose3(&mut r);
        }
        r
    }

    /// Rotation matrix as intervals; the interval mirror of [`Helmert7::rot_f64`].
    fn rot_iv(&self) -> [[Interval; 3]; 3] {
        let (f, t, p) = (pt(self.rx), pt(self.ry), pt(self.rz));
        let mut r = match self.rotation {
            Rotation::Exact => {
                let (cf, sf) = (f.cos(), f.sin());
                let (ct, st) = (t.cos(), t.sin());
                let (cp, sp) = (p.cos(), p.sin());
                [
                    [ct * cp, cf * sp + sf * st * cp, sf * sp - cf * st * cp],
                    [-ct * sp, cf * cp - sf * st * sp, sf * cp + cf * st * sp],
                    [st, -sf * ct, cf * ct],
                ]
            }
            Rotation::Approximate => [[pt(1.0), p, -t], [-p, pt(1.0), f], [t, -f, pt(1.0)]],
        };
        if self.convention == Convention::PositionVector {
            transpose3_iv(&mut r);
        }
        r
    }

    /// Ordinary `f64` 7-parameter Helmert forward: `out = s·R·X + T`, with
    /// `s = 1 + ds·10⁻⁶`. Matches PROJ `helmert_forward_3d` (`refp = 0`).
    pub fn forward_f64(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
        let r = self.rot_f64();
        let s = 1.0 + self.ds * 1e-6;
        let ox = s * (r[0][0] * x + r[0][1] * y + r[0][2] * z) + self.dx;
        let oy = s * (r[1][0] * x + r[1][1] * y + r[1][2] * z) + self.dy;
        let oz = s * (r[2][0] * x + r[2][1] * y + r[2][2] * z) + self.dz;
        (ox, oy, oz)
    }

    /// Guaranteed-bounds 7-parameter Helmert forward.
    ///
    /// Enclosure argument: the Helmert map is affine in the input (a scaled
    /// rotation plus translation). The rotation matrix is built once from the
    /// fixed parameters with enclosure primitives, and the matrix-vector product
    /// plus translation is a branch-free composition of `+`/`*`, so the result
    /// box provably encloses the exact shifted point.
    pub fn forward(&self, x: Interval, y: Interval, z: Interval) -> IntervalXyz {
        let r = self.rot_iv();
        let s = pt(1.0) + pt(self.ds) * pt(1e-6);
        let ox = s * (r[0][0] * x + r[0][1] * y + r[0][2] * z) + pt(self.dx);
        let oy = s * (r[1][0] * x + r[1][1] * y + r[1][2] * z) + pt(self.dy);
        let oz = s * (r[2][0] * x + r[2][1] * y + r[2][2] * z) + pt(self.dz);
        IntervalXyz::new(ox, oy, oz)
    }
}

#[inline]
fn transpose3(r: &mut [[f64; 3]; 3]) {
    let tmp = r[0][1];
    r[0][1] = r[1][0];
    r[1][0] = tmp;
    let tmp = r[0][2];
    r[0][2] = r[2][0];
    r[2][0] = tmp;
    let tmp = r[1][2];
    r[1][2] = r[2][1];
    r[2][1] = tmp;
}

#[inline]
fn transpose3_iv(r: &mut [[Interval; 3]; 3]) {
    let tmp = r[0][1];
    r[0][1] = r[1][0];
    r[1][0] = tmp;
    let tmp = r[0][2];
    r[0][2] = r[2][0];
    r[2][0] = tmp;
    let tmp = r[1][2];
    r[1][2] = r[2][1];
    r[2][1] = tmp;
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::const_params::{WGS84_A, WGS84_E, WGS84_ES};

    const DEG: f64 = core::f64::consts::PI / 180.0;

    fn merc_wgs84() -> MercParams {
        MercParams {
            a: WGS84_A,
            e: WGS84_E,
            k0: 1.0,
            x0: 0.0,
            y0: 0.0,
            lam0: 0.0,
        }
    }

    fn tmerc_wgs84(k0: f64) -> TmercParams {
        TmercParams {
            a: WGS84_A,
            es: WGS84_ES,
            k0,
            x0: 0.0,
            y0: 0.0,
            lam0: 0.0,
            phi0: 0.0,
        }
    }

    // ---- (a) point intervals contain the f64 result ----

    #[test]
    fn merc_point_interval_contains_f64() {
        let p = merc_wgs84();
        let samples = [
            (10.0, 45.0),
            (-73.5, 40.7),
            (0.0, 0.0),
            (179.0, -60.0),
            (5.0, 80.0),
        ];
        for &(lon, lat) in &samples {
            let (lam, phi) = (lon * DEG, lat * DEG);
            let (xf, yf) = p.forward_f64(lam, phi);
            let box_ = p.forward(Interval::point(lam), Interval::point(phi));
            assert!(
                box_.contains(xf, yf),
                "merc {lon},{lat}: {box_:?} !∋ ({xf},{yf})"
            );
        }
    }

    #[test]
    fn tmerc_point_interval_contains_f64() {
        let p = tmerc_wgs84(0.9996);
        let samples = [
            (3.0, 45.0),
            (1.0, 52.0),
            (0.0, 0.0),
            (2.5, -33.0),
            (-1.5, 70.0),
        ];
        for &(lon, lat) in &samples {
            let (lam, phi) = (lon * DEG, lat * DEG);
            let (xf, yf) = p.forward_f64(lam, phi);
            let box_ = p.forward(Interval::point(lam), Interval::point(phi));
            assert!(
                box_.contains(xf, yf),
                "tmerc {lon},{lat}: {box_:?} !∋ ({xf},{yf})"
            );
        }
    }

    #[test]
    fn helmert_point_interval_contains_f64() {
        for rotation in [Rotation::Approximate, Rotation::Exact] {
            for convention in [Convention::PositionVector, Convention::CoordinateFrame] {
                let h = Helmert7::from_arcsec_ppm(
                    10.0, 20.0, 30.0, 0.1, 0.2, 0.3, 1.0, convention, rotation,
                );
                let (x, y, z) = (4_000_000.0, 300_000.0, 5_000_000.0);
                let (xf, yf, zf) = h.forward_f64(x, y, z);
                let box_ = h.forward(Interval::point(x), Interval::point(y), Interval::point(z));
                assert!(
                    box_.contains(xf, yf, zf),
                    "helmert {rotation:?}/{convention:?}: {box_:?} !∋ ({xf},{yf},{zf})"
                );
            }
        }
    }

    #[test]
    fn point_interval_contains_f64_random_sweep() {
        // Deterministic pseudo-random sweep (LCG) for property (a).
        let merc = merc_wgs84();
        let tmerc = tmerc_wgs84(1.0);
        let mut state: u64 = 0x1234_5678_9abc_def0;
        let mut next = || {
            state = state
                .wrapping_mul(6364136223846793005)
                .wrapping_add(1442695040888963407);
            (state >> 11) as f64 / (1u64 << 53) as f64
        };
        for _ in 0..2000 {
            let lon = (next() * 160.0 - 80.0) * DEG;
            let lat = (next() * 140.0 - 70.0) * DEG;
            let (mxf, myf) = merc.forward_f64(lon, lat);
            let mb = merc.forward(Interval::point(lon), Interval::point(lat));
            assert!(mb.contains(mxf, myf), "merc sweep escape at {lon},{lat}");
            // tmerc: keep longitude near the central meridian (its valid domain).
            let tlon = (next() * 6.0 - 3.0) * DEG;
            let (txf, tyf) = tmerc.forward_f64(tlon, lat);
            let tb = tmerc.forward(Interval::point(tlon), Interval::point(lat));
            assert!(tb.contains(txf, tyf), "tmerc sweep escape at {tlon},{lat}");
        }
    }

    // ---- (b) monotone widening ----

    #[test]
    fn merc_monotone_widening() {
        let p = merc_wgs84();
        let lam0 = 10.0 * DEG;
        let phi0 = 45.0 * DEG;
        let narrow = p.forward(
            Interval::new(lam0 - 1e-4, lam0 + 1e-4),
            Interval::new(phi0 - 1e-4, phi0 + 1e-4),
        );
        let wide = p.forward(
            Interval::new(lam0 - 1e-2, lam0 + 1e-2),
            Interval::new(phi0 - 1e-2, phi0 + 1e-2),
        );
        assert!(wide.x.contains_interval(narrow.x), "merc x not monotone");
        assert!(wide.y.contains_interval(narrow.y), "merc y not monotone");
    }

    #[test]
    fn tmerc_monotone_widening() {
        let p = tmerc_wgs84(0.9996);
        let lam0 = 2.0 * DEG;
        let phi0 = 48.0 * DEG;
        let narrow = p.forward(
            Interval::new(lam0 - 1e-5, lam0 + 1e-5),
            Interval::new(phi0 - 1e-5, phi0 + 1e-5),
        );
        let wide = p.forward(
            Interval::new(lam0 - 1e-3, lam0 + 1e-3),
            Interval::new(phi0 - 1e-3, phi0 + 1e-3),
        );
        assert!(wide.x.contains_interval(narrow.x), "tmerc x not monotone");
        assert!(wide.y.contains_interval(narrow.y), "tmerc y not monotone");
    }

    #[test]
    fn helmert_monotone_widening() {
        let h = Helmert7::from_arcsec_ppm(
            10.0,
            20.0,
            30.0,
            0.1,
            0.2,
            0.3,
            1.0,
            Convention::PositionVector,
            Rotation::Exact,
        );
        let (x, y, z) = (4_000_000.0, 300_000.0, 5_000_000.0);
        let narrow = h.forward(
            Interval::new(x - 1.0, x + 1.0),
            Interval::new(y - 1.0, y + 1.0),
            Interval::new(z - 1.0, z + 1.0),
        );
        let wide = h.forward(
            Interval::new(x - 100.0, x + 100.0),
            Interval::new(y - 100.0, y + 100.0),
            Interval::new(z - 100.0, z + 100.0),
        );
        assert!(wide.x.contains_interval(narrow.x));
        assert!(wide.y.contains_interval(narrow.y));
        assert!(wide.z.contains_interval(narrow.z));
    }

    // ---- (c) known enclosures: agreement with reference PROJ 9.8.0 ----

    #[test]
    fn merc_matches_proj_reference() {
        // Empirical PROJ 9.8.0 (`proj +proj=merc +ellps=WGS84`, -f %.10f).
        let p = merc_wgs84();
        let cases = [
            (10.0, 45.0, 1_113_194.907_932_735, 5_591_295.918_553_391),
            (-73.5, 40.7, -8_181_982.573_305_607, 4_940_322.380_756_882),
        ];
        for &(lon, lat, ex, ey) in &cases {
            let (lam, phi) = (lon * DEG, lat * DEG);
            let (xf, yf) = p.forward_f64(lam, phi);
            assert!((xf - ex).abs() < 1e-4, "merc x {lon},{lat}: {xf} vs {ex}");
            assert!((yf - ey).abs() < 1e-4, "merc y {lon},{lat}: {yf} vs {ey}");
            // The certified box must also enclose the PROJ reference point.
            let box_ = p.forward(Interval::point(lam), Interval::point(phi));
            assert!(
                box_.contains(ex, ey),
                "merc box !∋ PROJ ref for {lon},{lat}"
            );
        }
    }

    #[test]
    fn tmerc_matches_proj_reference_submm() {
        // Empirical PROJ 9.8.0 (`proj +proj=tmerc +approx`). Our meridional arc
        // is the classic 6th-order-in-e series, so agreement is sub-millimetre;
        // we allow 5e-3 m (documented series-order tolerance).
        let cases = [
            (0.9996, 3.0, 45.0, 236_446.026_103_846, 4_987_329.504_699),
            (1.0, 1.0, 52.0, 68_677.175_765_481, 5_763_815.842_892_464),
        ];
        for &(k0, lon, lat, ex, ey) in &cases {
            let p = tmerc_wgs84(k0);
            let (lam, phi) = (lon * DEG, lat * DEG);
            let (xf, yf) = p.forward_f64(lam, phi);
            assert!((xf - ex).abs() < 5e-3, "tmerc x {lon},{lat}: {xf} vs {ex}");
            assert!((yf - ey).abs() < 5e-3, "tmerc y {lon},{lat}: {yf} vs {ey}");
        }
    }

    #[test]
    fn helmert_matches_proj_reference() {
        // Empirical PROJ 9.8.0 via `cct +proj=helmert ... -d 6`.
        let (x, y, z) = (4_000_000.0, 300_000.0, 5_000_000.0);
        // position_vector, approximate rotation.
        let pv = Helmert7::from_arcsec_ppm(
            10.0,
            20.0,
            30.0,
            0.1,
            0.2,
            0.3,
            1.0,
            Convention::PositionVector,
            Rotation::Approximate,
        );
        let (ox, oy, oz) = pv.forward_f64(x, y, z);
        assert!((ox - 4_000_018.411_809).abs() < 1e-4, "pv x {ox}");
        assert!((oy - 300_023.693_699).abs() < 1e-4, "pv y {oy}");
        assert!((oz - 5_000_031.266_931).abs() < 1e-4, "pv z {oz}");
        // coordinate_frame, approximate rotation.
        let cf = Helmert7::from_arcsec_ppm(
            10.0,
            20.0,
            30.0,
            0.1,
            0.2,
            0.3,
            1.0,
            Convention::CoordinateFrame,
            Rotation::Approximate,
        );
        let (cx, cy, cz) = cf.forward_f64(x, y, z);
        assert!((cx - 4_000_009.588_191).abs() < 1e-4, "cf x {cx}");
        assert!((cy - 300_016.906_301).abs() < 1e-4, "cf y {cy}");
        assert!((cz - 5_000_038.733_069).abs() < 1e-4, "cf z {cz}");
        // position_vector, exact rotation.
        let pe = Helmert7::from_arcsec_ppm(
            10.0,
            20.0,
            30.0,
            0.1,
            0.2,
            0.3,
            1.0,
            Convention::PositionVector,
            Rotation::Exact,
        );
        let (ex, ey, ez) = pe.forward_f64(x, y, z);
        assert!((ex - 4_000_018.411_803).abs() < 1e-4, "pe x {ex}");
        assert!((ey - 300_023.693_701).abs() < 1e-4, "pe y {ey}");
        assert!((ez - 5_000_031.266_931).abs() < 1e-4, "pe z {ez}");
    }

    #[test]
    fn known_rational_helmert_pure_translation() {
        // With zero rotation and zero scale the Helmert reduces to an exact
        // translation; the enclosure is tight and centred on the shifted point.
        let h = Helmert7::from_arcsec_ppm(
            100.0,
            -200.0,
            50.0,
            0.0,
            0.0,
            0.0,
            0.0,
            Convention::PositionVector,
            Rotation::Approximate,
        );
        let box_ = h.forward(
            Interval::point(1000.0),
            Interval::point(2000.0),
            Interval::point(3000.0),
        );
        assert!(box_.contains(1100.0, 1800.0, 3050.0));
        assert!(box_.x.width() < 1e-6 && box_.y.width() < 1e-6 && box_.z.width() < 1e-6);
    }

    #[test]
    fn merc_equator_is_origin() {
        // Known enclosure: the equator maps to y = 0 exactly (isometric lat 0).
        let p = merc_wgs84();
        let box_ = p.forward(Interval::point(0.0), Interval::point(0.0));
        assert!(box_.contains(0.0, 0.0));
        assert!(box_.y.width() < 1e-6);
    }

    // ---- (d) adversarial enclosure: the box from a WIDE input interval must
    // contain the exact f64 forward at EVERY interior sample, not merely the
    // centre. This is the core T7.4 soundness contract; a too-tight interval
    // (broken outward rounding, or a non-monotone step evaluated only at the
    // endpoints) would let an interior point escape the box. ----

    #[test]
    fn merc_box_encloses_all_interior_points() {
        let p = merc_wgs84();
        let lam_iv = Interval::new(5.0 * DEG, 15.0 * DEG);
        let phi_iv = Interval::new(30.0 * DEG, 60.0 * DEG);
        let box_ = p.forward(lam_iv, phi_iv);
        let (nl, np) = (37usize, 41usize);
        for i in 0..=nl {
            let lam = lam_iv.lo() + (lam_iv.hi() - lam_iv.lo()) * i as f64 / nl as f64;
            for j in 0..=np {
                let phi = phi_iv.lo() + (phi_iv.hi() - phi_iv.lo()) * j as f64 / np as f64;
                let (xf, yf) = p.forward_f64(lam, phi);
                assert!(
                    box_.contains(xf, yf),
                    "merc interior ({lam},{phi}) escaped {box_:?}: ({xf},{yf})"
                );
            }
        }
    }

    #[test]
    fn tmerc_box_encloses_all_interior_points() {
        let p = tmerc_wgs84(0.9996);
        // Keep longitude within the transverse-Mercator valid band around the
        // central meridian; sweep a wide latitude band across the equator.
        let lam_iv = Interval::new(-2.5 * DEG, 2.5 * DEG);
        let phi_iv = Interval::new(-40.0 * DEG, 55.0 * DEG);
        let box_ = p.forward(lam_iv, phi_iv);
        let (nl, np) = (33usize, 61usize);
        for i in 0..=nl {
            let lam = lam_iv.lo() + (lam_iv.hi() - lam_iv.lo()) * i as f64 / nl as f64;
            for j in 0..=np {
                let phi = phi_iv.lo() + (phi_iv.hi() - phi_iv.lo()) * j as f64 / np as f64;
                let (xf, yf) = p.forward_f64(lam, phi);
                assert!(
                    box_.contains(xf, yf),
                    "tmerc interior ({lam},{phi}) escaped {box_:?}: ({xf},{yf})"
                );
            }
        }
    }

    #[test]
    fn helmert_box_encloses_all_interior_points() {
        for rotation in [Rotation::Approximate, Rotation::Exact] {
            for convention in [Convention::PositionVector, Convention::CoordinateFrame] {
                let h = Helmert7::from_arcsec_ppm(
                    10.0, 20.0, 30.0, 0.4, -0.6, 0.9, 3.0, convention, rotation,
                );
                let x_iv = Interval::new(3_900_000.0, 4_100_000.0);
                let y_iv = Interval::new(250_000.0, 350_000.0);
                let z_iv = Interval::new(4_900_000.0, 5_100_000.0);
                let box_ = h.forward(x_iv, y_iv, z_iv);
                let n = 7usize;
                for i in 0..=n {
                    let x = x_iv.lo() + (x_iv.hi() - x_iv.lo()) * i as f64 / n as f64;
                    for j in 0..=n {
                        let y = y_iv.lo() + (y_iv.hi() - y_iv.lo()) * j as f64 / n as f64;
                        for k in 0..=n {
                            let z = z_iv.lo() + (z_iv.hi() - z_iv.lo()) * k as f64 / n as f64;
                            let (xf, yf, zf) = h.forward_f64(x, y, z);
                            assert!(
                                box_.contains(xf, yf, zf),
                                "helmert {rotation:?}/{convention:?} interior \
                                 ({x},{y},{z}) escaped {box_:?}"
                            );
                        }
                    }
                }
            }
        }
    }
}