oxiproj-transformations 0.1.2

//! `helmert` and `molobadekas` datum transformations.
//!
//! Faithful Pure-Rust port of PROJ `src/transformations/helmert.cpp`.
//!
//! Implements 3(6)-, 4(8)- and 7(14)-parameter Helmert transformations for 3D
//! data, plus the Molodensky-Badekas variant where the rotation is applied
//! about a reference point (pivot) rather than the centre of the spheroid.
//!
//! Translations are in metres, rotations in arc-seconds (converted to radians),
//! scale in parts-per-million. The optional `t_epoch`/rate parameters give the
//! kinematic (time-dependent) counterparts.
//!
//! `towgs84` shorthand parsing and the (now-rejected) `transpose` flag are
//! both implemented below, matching PROJ 9.8.0's `read_convention` and
//! `init_helmert_six_parameters`. The one remaining deferred item (not yet
//! ported from `helmert.cpp`) is the legacy `+axis neu` handling.

use crate::{TransBuild, TransParamLookup, TransParams};
use oxiproj_core::{Coord, IoUnits, Lpz, Operation, ProjError, ProjResult, Xyz, DEG_TO_RAD};

/// Arc-seconds to radians. Mirrors PROJ's per-arcsecond scaling in
/// `src/transformations/helmert.cpp`.
const ARCSEC_TO_RAD: f64 = DEG_TO_RAD / 3600.0;

/// Opaque Helmert state, mirroring `pj_opaque_helmert` in
/// `src/transformations/helmert.cpp`.
#[derive(Debug)]
struct Helmert {
    /// Base translation (metres).
    xyz_0: [f64; 3],
    /// Translation rates (metres/year).
    dxyz: [f64; 3],
    /// Molodensky-Badekas pivot (metres); `[0, 0, 0]` for plain Helmert.
    refp: [f64; 3],
    /// Base rotations (omega, phi, kappa) in radians.
    opk_0: [f64; 3],
    /// Rotation rates (radians/year).
    dopk: [f64; 3],
    /// Base scale (ppm).
    scale_0: f64,
    /// Scale rate (ppm/year).
    dscale: f64,
    /// Base rotation angle for the 2D 4-parameter case (radians).
    theta_0: f64,
    /// Rotation-angle rate for the 2D 4-parameter case (radians/year).
    dtheta: f64,
    /// Central epoch for the kinematic parameters.
    t_epoch: f64,
    /// True when no rotation is requested (rotations and their rates all zero).
    no_rotation: bool,
    /// True to use the exact (non-linearised) rotation matrix.
    exact: bool,
    /// True for the 2D 4-parameter (theta/scale) variant.
    fourparam: bool,
    /// True for the position-vector convention, false for coordinate-frame.
    is_position_vector: bool,
}

/// Epoch-derived parameters, mirroring the locals updated by
/// `update_parameters` in `src/transformations/helmert.cpp`.
struct Derived {
    /// Epoch-adjusted translation (metres).
    xyz: [f64; 3],
    /// Epoch-adjusted rotations (radians).
    opk: [f64; 3],
    /// Epoch-adjusted scale (ppm).
    scale: f64,
    /// Epoch-adjusted 2D rotation angle (radians).
    theta: f64,
}

/// Build the rotation matrix, ported from `build_rot_matrix` in
/// `src/transformations/helmert.cpp`.
///
/// The equations are given in the coordinate-frame convention; passing
/// `is_pv = true` transposes the off-diagonal terms for the position-vector
/// convention.
fn rot_matrix(opk: [f64; 3], exact: bool, is_pv: bool) -> [[f64; 3]; 3] {
    let f = opk[0];
    let t = opk[1];
    let p = opk[2];

    let mut r = if exact {
        let cf = f.cos();
        let sf = f.sin();
        let ct = t.cos();
        let st = t.sin();
        let cp = p.cos();
        let sp = 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],
        ]
    } else {
        [[1.0, p, -t], [-p, 1.0, f], [t, -f, 1.0]]
    };

    if is_pv {
        let tmp01 = r[0][1];
        r[0][1] = r[1][0];
        r[1][0] = tmp01;
        let tmp02 = r[0][2];
        r[0][2] = r[2][0];
        r[2][0] = tmp02;
        let tmp12 = r[1][2];
        r[1][2] = r[2][1];
        r[2][1] = tmp12;
    }
    r
}

impl Helmert {
    /// Update the base parameters to the observation epoch `t_obs`.
    /// Ported from `update_parameters` in `src/transformations/helmert.cpp`.
    fn derive(&self, t_obs: f64) -> Derived {
        let dt = t_obs - self.t_epoch;
        let xyz = core::array::from_fn(|i| self.xyz_0[i] + self.dxyz[i] * dt);
        let opk = core::array::from_fn(|i| self.opk_0[i] + self.dopk[i] * dt);
        Derived {
            xyz,
            opk,
            scale: self.scale_0 + self.dscale * dt,
            theta: self.theta_0 + self.dtheta * dt,
        }
    }

    /// Forward 3D core. Ported from `helmert_forward_3d` in
    /// `src/transformations/helmert.cpp`. `x`, `y`, `z` are the input
    /// Cartesian coordinates.
    // ported verbatim from PROJ src/transformations/helmert.cpp
    #[allow(clippy::float_cmp)]
    fn fwd_core(&self, x: f64, y: f64, z: f64, d: &Derived, r: &[[f64; 3]; 3]) -> (f64, f64, f64) {
        if self.fourparam {
            let cr = d.theta.cos() * d.scale;
            let sr = d.theta.sin() * d.scale;
            let ox = cr * x + sr * y + self.xyz_0[0];
            let oy = -sr * x + cr * y + self.xyz_0[1];
            return (ox, oy, z);
        }
        if self.no_rotation && d.scale == 0.0 {
            return (x + d.xyz[0], y + d.xyz[1], z + d.xyz[2]);
        }
        let scale_factor = 1.0 + d.scale * 1e-6;
        let xr = x - self.refp[0];
        let yr = y - self.refp[1];
        let zr = z - self.refp[2];
        let ox = scale_factor * (r[0][0] * xr + r[0][1] * yr + r[0][2] * zr) + d.xyz[0];
        let oy = scale_factor * (r[1][0] * xr + r[1][1] * yr + r[1][2] * zr) + d.xyz[1];
        let oz = scale_factor * (r[2][0] * xr + r[2][1] * yr + r[2][2] * zr) + d.xyz[2];
        (ox, oy, oz)
    }

    /// Reverse 3D core. Ported from `helmert_reverse_3d` in
    /// `src/transformations/helmert.cpp`. `x`, `y`, `z` are the input
    /// Cartesian coordinates (the forward output).
    // ported verbatim from PROJ src/transformations/helmert.cpp
    #[allow(clippy::float_cmp)]
    fn inv_core(&self, x: f64, y: f64, z: f64, d: &Derived, r: &[[f64; 3]; 3]) -> (f64, f64, f64) {
        if self.fourparam {
            let cr = d.theta.cos() / d.scale;
            let sr = d.theta.sin() / d.scale;
            let xx = x - self.xyz_0[0];
            let yy = y - self.xyz_0[1];
            let ox = xx * cr - yy * sr;
            let oy = xx * sr + yy * cr;
            return (ox, oy, z);
        }
        if self.no_rotation && d.scale == 0.0 {
            return (x - d.xyz[0], y - d.xyz[1], z - d.xyz[2]);
        }
        let scale_factor = 1.0 + d.scale * 1e-6;
        let xr = (x - d.xyz[0]) / scale_factor;
        let yr = (y - d.xyz[1]) / scale_factor;
        let zr = (z - d.xyz[2]) / scale_factor;
        let ox = (r[0][0] * xr + r[1][0] * yr + r[2][0] * zr) + self.refp[0];
        let oy = (r[0][1] * xr + r[1][1] * yr + r[2][1] * zr) + self.refp[1];
        let oz = (r[0][2] * xr + r[1][2] * yr + r[2][2] * zr) + self.refp[2];
        (ox, oy, oz)
    }
}

impl Operation for Helmert {
    fn forward_3d(&self, lpz: Lpz) -> ProjResult<Xyz> {
        let d = self.derive(self.t_epoch);
        let r = rot_matrix(d.opk, self.exact, self.is_position_vector);
        let (ox, oy, oz) = self.fwd_core(lpz.lam, lpz.phi, lpz.z, &d, &r);
        Ok(Xyz::new(ox, oy, oz))
    }

    fn inverse_3d(&self, xyz: Xyz) -> ProjResult<Lpz> {
        let d = self.derive(self.t_epoch);
        let r = rot_matrix(d.opk, self.exact, self.is_position_vector);
        let (ox, oy, oz) = self.inv_core(xyz.x, xyz.y, xyz.z, &d, &r);
        Ok(Lpz::new(ox, oy, oz))
    }

    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let t = v[3];
        let t_obs = if t.is_finite() { t } else { self.t_epoch };
        let d = self.derive(t_obs);
        let r = rot_matrix(d.opk, self.exact, self.is_position_vector);
        let (ox, oy, oz) = self.fwd_core(v[0], v[1], v[2], &d, &r);
        Ok(Coord::new(ox, oy, oz, t))
    }

    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let t = v[3];
        let t_obs = if t.is_finite() { t } else { self.t_epoch };
        let d = self.derive(t_obs);
        let r = rot_matrix(d.opk, self.exact, self.is_position_vector);
        let (ox, oy, oz) = self.inv_core(v[0], v[1], v[2], &d, &r);
        Ok(Coord::new(ox, oy, oz, t))
    }

    fn has_inverse(&self) -> bool {
        true
    }
}

/// Resolve the rotation convention. Ported from `read_convention` in
/// `src/transformations/helmert.cpp`: when a rotation is present the
/// `convention` keyword is mandatory — `+towgs84` does **not** exempt the
/// caller from supplying it, but it does restrict the accepted value.
///
/// Empirically verified against `cct` (PROJ 9.7.0, same `read_convention`
/// logic in 9.8.0):
/// * `+proj=helmert +towgs84=1,2,3,0.1,0.2,0.3,0.5` (no `+convention`) →
///   error 1026 "missing 'convention' argument".
/// * `... +convention=coordinate_frame` → error 1027 "towgs84 should only be
///   used with convention=position_vector".
/// * `... +convention=position_vector` → succeeds.
fn read_convention(
    params: &dyn TransParamLookup,
    no_rotation: bool,
    towgs84_present: bool,
) -> ProjResult<bool> {
    if no_rotation {
        return Ok(false);
    }
    let is_position_vector = match params.get_str("convention") {
        Some("position_vector") => true,
        Some("coordinate_frame") => false,
        None => return Err(ProjError::MissingArg),
        Some(_) => return Err(ProjError::IllegalArgValue),
    };
    // PROJ: "historically towgs84 in PROJ has always been using
    // position_vector convention. Accepting coordinate_frame would be
    // confusing."
    if towgs84_present && !is_position_vector {
        return Err(ProjError::IllegalArgValue);
    }
    Ok(is_position_vector)
}

/// Construct the `helmert` transformation. Ported from
/// `PJ_TRANSFORMATION(helmert, ...)` in `src/transformations/helmert.cpp`.
// ported verbatim from PROJ src/transformations/helmert.cpp
#[allow(clippy::float_cmp)]
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let params = p.params;

    let mut xyz_0 = [
        params.get_f64("x").unwrap_or(0.0),
        params.get_f64("y").unwrap_or(0.0),
        params.get_f64("z").unwrap_or(0.0),
    ];
    let mut opk_0 = [
        params.get_f64("rx").unwrap_or(0.0) * ARCSEC_TO_RAD,
        params.get_f64("ry").unwrap_or(0.0) * ARCSEC_TO_RAD,
        params.get_f64("rz").unwrap_or(0.0) * ARCSEC_TO_RAD,
    ];
    let exact = params.get_bool("exact");

    let has_theta = params.exists("theta");
    let (left, right) = if has_theta {
        (IoUnits::Projected, IoUnits::Projected)
    } else {
        (IoUnits::Cartesian, IoUnits::Cartesian)
    };

    let mut fourparam = false;
    let mut theta_0 = 0.0;
    let mut scale_0: f64 = 0.0;
    if has_theta {
        theta_0 = params.get_f64("theta").unwrap_or(0.0) * ARCSEC_TO_RAD;
        fourparam = true;
        scale_0 = 1.0;
    }
    if params.exists("s") {
        scale_0 = params.get_f64("s").unwrap_or(0.0);
        if scale_0 <= -1.0e6 {
            return Err(ProjError::IllegalArgValue);
        }
        if fourparam && scale_0 == 0.0 {
            return Err(ProjError::IllegalArgValue);
        }
    }

    let dxyz = [
        params.get_f64("dx").unwrap_or(0.0),
        params.get_f64("dy").unwrap_or(0.0),
        params.get_f64("dz").unwrap_or(0.0),
    ];
    let dopk = [
        params.get_f64("drx").unwrap_or(0.0) * ARCSEC_TO_RAD,
        params.get_f64("dry").unwrap_or(0.0) * ARCSEC_TO_RAD,
        params.get_f64("drz").unwrap_or(0.0) * ARCSEC_TO_RAD,
    ];
    let dtheta = params.get_f64("dtheta").unwrap_or(0.0) * ARCSEC_TO_RAD;
    let dscale = params.get_f64("ds").unwrap_or(0.0);
    let t_epoch = params.get_f64("t_epoch").unwrap_or(0.0);

    // Detect deprecated `transpose` flag and return an error.
    // Matches PROJ 9.8.0 behaviour: the flag was removed in favour of
    // the explicit `convention=position_vector/coordinate_frame` parameter.
    if params.get_bool("transpose") {
        return Err(ProjError::IllegalArgValue);
    }

    // Support the classic PROJ `+towgs84=dx,dy,dz[,rx,ry,rz,s]` shorthand.
    // In PROJ's C code this is parsed at init time and stored in datum_params[];
    // here we parse it directly from the string parameter.
    let towgs84_present = params.exists("towgs84");
    if let Some(tw) = params.get_str("towgs84") {
        let vals: Vec<f64> = tw
            .split(',')
            .map(|s| s.trim().parse::<f64>())
            .collect::<Result<Vec<_>, _>>()
            .map_err(|_| ProjError::IllegalArgValue)?;
        match vals.len() {
            3 => {
                // 3-parameter shift: dx, dy, dz
                xyz_0[0] = vals[0];
                xyz_0[1] = vals[1];
                xyz_0[2] = vals[2];
            }
            7 => {
                // 7-parameter shift: dx, dy, dz, rx, ry, rz, s
                xyz_0[0] = vals[0];
                xyz_0[1] = vals[1];
                xyz_0[2] = vals[2];
                opk_0[0] = vals[3] * ARCSEC_TO_RAD;
                opk_0[1] = vals[4] * ARCSEC_TO_RAD;
                opk_0[2] = vals[5] * ARCSEC_TO_RAD;
                // PROJ historically stores scale as ppm; towgs84 scale IS ppm
                scale_0 = vals[6];
            }
            _ => return Err(ProjError::IllegalArgValue),
        }
    }

    // `no_rotation` must be computed after towgs84 may have updated opk_0.
    let no_rotation = opk_0.iter().all(|&v| v == 0.0) && dopk.iter().all(|&v| v == 0.0);
    let is_position_vector = read_convention(params, no_rotation, towgs84_present)?;

    let op = Helmert {
        xyz_0,
        dxyz,
        refp: [0.0, 0.0, 0.0],
        opk_0,
        dopk,
        scale_0,
        dscale,
        theta_0,
        dtheta,
        t_epoch,
        no_rotation,
        exact,
        fourparam,
        is_position_vector,
    };
    Ok(TransBuild::new(Box::new(op), left, right))
}

/// Construct the `molobadekas` transformation. Ported from
/// `PJ_TRANSFORMATION(molobadekas, ...)` in `src/transformations/helmert.cpp`.
// ported verbatim from PROJ src/transformations/helmert.cpp
#[allow(clippy::float_cmp)]
pub fn new_molobadekas(p: &TransParams) -> ProjResult<TransBuild> {
    let params = p.params;

    let mut xyz_0 = [
        params.get_f64("x").unwrap_or(0.0),
        params.get_f64("y").unwrap_or(0.0),
        params.get_f64("z").unwrap_or(0.0),
    ];
    let opk_0 = [
        params.get_f64("rx").unwrap_or(0.0) * ARCSEC_TO_RAD,
        params.get_f64("ry").unwrap_or(0.0) * ARCSEC_TO_RAD,
        params.get_f64("rz").unwrap_or(0.0) * ARCSEC_TO_RAD,
    ];
    let exact = params.get_bool("exact");
    let scale_0 = params.get_f64("s").unwrap_or(0.0);

    let dxyz = [0.0, 0.0, 0.0];
    let dopk = [0.0, 0.0, 0.0];

    // PROJ's `PJ_TRANSFORMATION(molobadekas, ...)` leaves `no_rotation` at its
    // calloc'd value of 0 (false) and never runs the rotation-is-zero check used
    // by `helmert`. The Molodensky-Badekas forward/reverse therefore always take
    // the full `scale * R * (x - refp) + xyz_0` path so the pivot offset is
    // honoured even when the rotation happens to be zero. Consequently the
    // `convention` keyword is always mandatory here.
    let no_rotation = false;
    // Molodensky-Badekas has no `+towgs84` shorthand in PROJ (only the plain
    // `helmert` operator's init function parses it), so the towgs84/
    // convention restriction never applies here.
    let is_position_vector = read_convention(params, no_rotation, false)?;

    let refp = [
        params.get_f64("px").unwrap_or(0.0),
        params.get_f64("py").unwrap_or(0.0),
        params.get_f64("pz").unwrap_or(0.0),
    ];
    // PROJ optimisation: fold the pivot offset into the translation so the
    // forward core can subtract `refp` before rotation and add it back via the
    // translation term.
    for i in 0..3 {
        xyz_0[i] += refp[i];
    }

    let op = Helmert {
        xyz_0,
        dxyz,
        refp,
        opk_0,
        dopk,
        scale_0,
        dscale: 0.0,
        theta_0: 0.0,
        dtheta: 0.0,
        t_epoch: 0.0,
        no_rotation,
        exact,
        fourparam: false,
        is_position_vector,
    };
    Ok(TransBuild::new(
        Box::new(op),
        IoUnits::Cartesian,
        IoUnits::Cartesian,
    ))
}

#[cfg(test)]
mod tests {
    use super::*;
    use oxiproj_core::Ellipsoid;
    use std::collections::HashMap;

    struct Lookup {
        nums: HashMap<&'static str, f64>,
        convention: Option<&'static str>,
        exact: bool,
    }

    impl Lookup {
        fn new(pairs: &[(&'static str, f64)]) -> Lookup {
            let mut nums = HashMap::new();
            for &(k, v) in pairs {
                nums.insert(k, v);
            }
            Lookup {
                nums,
                convention: None,
                exact: false,
            }
        }
        fn with_convention(mut self, c: &'static str) -> Lookup {
            self.convention = Some(c);
            self
        }
        fn with_exact(mut self, e: bool) -> Lookup {
            self.exact = e;
            self
        }
    }

    impl TransParamLookup for Lookup {
        fn get_dms(&self, _key: &str) -> Option<f64> {
            None
        }
        fn get_f64(&self, key: &str) -> Option<f64> {
            self.nums.get(key).copied()
        }
        fn get_int(&self, _key: &str) -> Option<i64> {
            None
        }
        fn get_str(&self, key: &str) -> Option<&str> {
            if key == "convention" {
                self.convention
            } else {
                None
            }
        }
        fn get_bool(&self, key: &str) -> bool {
            key == "exact" && self.exact
        }
        fn exists(&self, key: &str) -> bool {
            if key == "convention" {
                self.convention.is_some()
            } else {
                self.nums.contains_key(key)
            }
        }
    }

    fn wgs84() -> Ellipsoid {
        Ellipsoid::named("WGS84").unwrap()
    }

    fn close(a: f64, b: f64, tol: f64, label: &str) {
        assert!(
            (a - b).abs() < tol,
            "{label}: {a} vs {b} (diff {})",
            (a - b).abs()
        );
    }

    #[test]
    fn coordinate_frame_seven_param() {
        let ell = wgs84();
        let lk = Lookup::new(&[
            ("x", 0.06155),
            ("y", -0.01087),
            ("z", -0.04019),
            ("rx", -0.0394924),
            ("ry", -0.0327221),
            ("rz", -0.0328979),
            ("s", -0.009994),
        ])
        .with_convention("coordinate_frame");
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk,
            registry: None,
        };
        let b = new(&pp).unwrap();

        let out = b
            .operation
            .forward_3d(Lpz::new(3565285.0, 855949.0, 5201383.0))
            .unwrap();
        close(out.x, 3565285.714553948, 1e-6, "cf.x");
        close(out.y, 855948.553335006, 1e-6, "cf.y");
        close(out.z, 5201382.506110126, 1e-6, "cf.z");

        let out4 = b
            .operation
            .forward_4d(Coord::new(3565285.0, 855949.0, 5201383.0, 0.0))
            .unwrap();
        let v = out4.v();
        close(v[0], out.x, 1e-9, "cf4.x");
        close(v[1], out.y, 1e-9, "cf4.y");
        close(v[2], out.z, 1e-9, "cf4.z");

        let back = b.operation.inverse_3d(out).unwrap();
        close(back.lam, 3565285.0, 1e-6, "cf.inv.x");
        close(back.phi, 855949.0, 1e-6, "cf.inv.y");
        close(back.z, 5201383.0, 1e-6, "cf.inv.z");
    }

    #[test]
    fn position_vector_seven_param() {
        let ell = wgs84();
        let lk = Lookup::new(&[
            ("x", 0.06155),
            ("y", -0.01087),
            ("z", -0.04019),
            ("rx", -0.0394924),
            ("ry", -0.0327221),
            ("rz", -0.0328979),
            ("s", -0.009994),
        ])
        .with_convention("position_vector");
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk,
            registry: None,
        };
        let b = new(&pp).unwrap();

        let out = b
            .operation
            .forward_3d(Lpz::new(3565285.0, 855949.0, 5201383.0))
            .unwrap();
        close(out.x, 3565284.337283136, 1e-6, "pv.x");
        close(out.y, 855949.407816285, 1e-6, "pv.y");
        close(out.z, 5201383.30954463, 1e-6, "pv.z");

        let back = b.operation.inverse_3d(out).unwrap();
        close(back.lam, 3565285.0, 1e-6, "pv.inv.x");
        close(back.phi, 855949.0, 1e-6, "pv.inv.y");
        close(back.z, 5201383.0, 1e-6, "pv.inv.z");
    }

    #[test]
    fn coordinate_frame_exact() {
        let ell = wgs84();
        let lk = Lookup::new(&[
            ("x", 0.06155),
            ("y", -0.01087),
            ("z", -0.04019),
            ("rx", -0.0394924),
            ("ry", -0.0327221),
            ("rz", -0.0328979),
            ("s", -0.009994),
        ])
        .with_convention("coordinate_frame")
        .with_exact(true);
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk,
            registry: None,
        };
        let b = new(&pp).unwrap();

        let out = b
            .operation
            .forward_3d(Lpz::new(3565285.0, 855949.0, 5201383.0))
            .unwrap();
        close(out.x, 3565285.714554043, 1e-6, "ex.x");
        close(out.y, 855948.553335111, 1e-6, "ex.y");
        close(out.z, 5201382.506109965, 1e-6, "ex.z");

        let back = b.operation.inverse_3d(out).unwrap();
        close(back.lam, 3565285.0, 1e-6, "ex.inv.x");
        close(back.phi, 855949.0, 1e-6, "ex.inv.y");
        close(back.z, 5201383.0, 1e-6, "ex.inv.z");
    }

    #[test]
    fn three_param_translation() {
        let ell = wgs84();
        let lk = Lookup::new(&[("x", 10.0), ("y", 20.0), ("z", 30.0)]);
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk,
            registry: None,
        };
        let b = new(&pp).unwrap();

        let out = b
            .operation
            .forward_3d(Lpz::new(3565285.0, 855949.0, 5201383.0))
            .unwrap();
        close(out.x, 3565295.0, 1e-9, "t3.x");
        close(out.y, 855969.0, 1e-9, "t3.y");
        close(out.z, 5201413.0, 1e-9, "t3.z");

        let back = b.operation.inverse_3d(out).unwrap();
        close(back.lam, 3565285.0, 1e-9, "t3.inv.x");
        close(back.phi, 855949.0, 1e-9, "t3.inv.y");
        close(back.z, 5201383.0, 1e-9, "t3.inv.z");
    }

    #[test]
    fn time_dependent_fifteen_param() {
        let ell = wgs84();
        let lk = Lookup::new(&[
            ("x", 0.0521),
            ("y", 0.0493),
            ("z", -0.0585),
            ("rx", 0.000891),
            ("ry", 0.00539),
            ("rz", -0.008712),
            ("s", 0.001034),
            ("dx", 0.0001),
            ("dy", 0.0001),
            ("dz", -0.0018),
            ("drx", 0.000081),
            ("dry", 0.000490),
            ("drz", -0.000792),
            ("ds", 0.000108),
            ("t_epoch", 2000.0),
        ])
        .with_convention("position_vector");
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk,
            registry: None,
        };
        let b = new(&pp).unwrap();

        let out2010 = b
            .operation
            .forward_4d(Coord::new(3565285.0, 855949.0, 5201383.0, 2010.0))
            .unwrap();
        let v = out2010.v();
        close(v[0], 3565285.389138885, 1e-6, "td.x");
        close(v[1], 855948.721731386, 1e-6, "td.y");
        close(v[2], 5201382.76369192, 1e-6, "td.z");

        let out2000 = b
            .operation
            .forward_4d(Coord::new(3565285.0, 855949.0, 5201383.0, 2000.0))
            .unwrap();
        let w = out2000.v();
        let dx = (v[0] - w[0]).abs();
        let dy = (v[1] - w[1]).abs();
        assert!(
            dx > 1e-4 || dy > 1e-4,
            "time must change result: dx={dx} dy={dy}"
        );

        let back = b.operation.inverse_4d(out2010).unwrap();
        let bk = back.v();
        close(bk[0], 3565285.0, 1e-6, "td.inv.x");
        close(bk[1], 855949.0, 1e-6, "td.inv.y");
        close(bk[2], 5201383.0, 1e-6, "td.inv.z");
    }

    #[test]
    fn molobadekas_identity() {
        let ell = wgs84();
        let lk = Lookup::new(&[
            ("x", 0.0),
            ("y", 0.0),
            ("z", 0.0),
            ("rx", 0.0),
            ("ry", 0.0),
            ("rz", 0.0),
            ("s", 0.0),
            ("px", 3789470.71),
            ("py", -115924.04),
            ("pz", 5111432.06),
        ])
        .with_convention("coordinate_frame");
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk,
            registry: None,
        };
        let b = new_molobadekas(&pp).unwrap();

        let out = b
            .operation
            .forward_3d(Lpz::new(2845456.0813, 2160954.2453, 5265993.2296))
            .unwrap();
        close(out.x, 2845456.0813, 1e-6, "mb.x");
        close(out.y, 2160954.2453, 1e-6, "mb.y");
        close(out.z, 5265993.2296, 1e-6, "mb.z");

        let back = b.operation.inverse_3d(out).unwrap();
        close(back.lam, 2845456.0813, 1e-6, "mb.inv.x");
        close(back.phi, 2160954.2453, 1e-6, "mb.inv.y");
        close(back.z, 5265993.2296, 1e-6, "mb.inv.z");
    }

    #[test]
    fn molobadekas_pivot_round_trip() {
        let ell = wgs84();
        // Use the exact rotation matrix so its transpose is a true inverse and
        // the forward/reverse pair round-trips to near machine precision; the
        // linearised matrix only reaches ~1e-5 m at these coordinate magnitudes.
        let lk = Lookup::new(&[
            ("rx", 1.0),
            ("ry", -2.0),
            ("rz", 0.5),
            ("s", 2.0),
            ("px", 3789470.71),
            ("py", -115924.04),
            ("pz", 5111432.06),
        ])
        .with_convention("coordinate_frame")
        .with_exact(true);
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk,
            registry: None,
        };
        let b = new_molobadekas(&pp).unwrap();

        let input = Lpz::new(2845456.0813, 2160954.2453, 5265993.2296);
        let out = b.operation.forward_3d(input).unwrap();
        let back = b
            .operation
            .inverse_3d(Xyz::new(out.x, out.y, out.z))
            .unwrap();
        close(back.lam, input.lam, 1e-6, "mbp.inv.x");
        close(back.phi, input.phi, 1e-6, "mbp.inv.y");
        close(back.z, input.z, 1e-6, "mbp.inv.z");
    }

    #[test]
    fn convention_errors() {
        let ell = wgs84();

        let lk_missing = Lookup::new(&[("rx", 1.0), ("ry", 2.0), ("rz", 3.0)]);
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk_missing,
            registry: None,
        };
        assert_eq!(new(&pp).err(), Some(ProjError::MissingArg));

        let lk_bad = Lookup::new(&[("rx", 1.0), ("ry", 2.0), ("rz", 3.0)]).with_convention("foo");
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk_bad,
            registry: None,
        };
        assert_eq!(new(&pp).err(), Some(ProjError::IllegalArgValue));
    }

    /// Test helper: supplies a `towgs84` string parameter and, optionally, an
    /// explicit `convention`.
    struct LookupWithTowgs84 {
        towgs84: String,
        convention: Option<&'static str>,
    }

    impl LookupWithTowgs84 {
        fn new(val: &str) -> Self {
            Self {
                towgs84: val.to_owned(),
                convention: None,
            }
        }
        fn with_convention(mut self, c: &'static str) -> Self {
            self.convention = Some(c);
            self
        }
    }

    impl TransParamLookup for LookupWithTowgs84 {
        fn get_dms(&self, _key: &str) -> Option<f64> {
            None
        }
        fn get_f64(&self, _key: &str) -> Option<f64> {
            None
        }
        fn get_int(&self, _key: &str) -> Option<i64> {
            None
        }
        fn get_str(&self, key: &str) -> Option<&str> {
            match key {
                "towgs84" => Some(&self.towgs84),
                "convention" => self.convention,
                _ => None,
            }
        }
        fn get_bool(&self, _key: &str) -> bool {
            false
        }
        fn exists(&self, key: &str) -> bool {
            match key {
                "towgs84" => true,
                "convention" => self.convention.is_some(),
                _ => false,
            }
        }
    }

    /// Test helper: pretends that `+transpose` is present.
    struct LookupWithTranspose;

    impl LookupWithTranspose {
        fn new() -> Self {
            Self
        }
    }

    impl TransParamLookup for LookupWithTranspose {
        fn get_dms(&self, _key: &str) -> Option<f64> {
            None
        }
        fn get_f64(&self, _key: &str) -> Option<f64> {
            None
        }
        fn get_int(&self, _key: &str) -> Option<i64> {
            None
        }
        fn get_str(&self, _key: &str) -> Option<&str> {
            None
        }
        fn get_bool(&self, key: &str) -> bool {
            key == "transpose"
        }
        fn exists(&self, _key: &str) -> bool {
            false
        }
    }

    #[test]
    fn towgs84_three_param() {
        let ell = wgs84();
        // +towgs84=1.0,2.0,3.0 should shift XYZ by (1, 2, 3)
        let lk = LookupWithTowgs84::new("1.0,2.0,3.0");
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk,
            registry: None,
        };
        let b = new(&pp).unwrap();
        let out = b
            .operation
            .forward_3d(Lpz::new(100.0, 200.0, 300.0))
            .unwrap();
        close(out.x, 101.0, 1e-9, "tw3.x");
        close(out.y, 202.0, 1e-9, "tw3.y");
        close(out.z, 303.0, 1e-9, "tw3.z");

        let back = b.operation.inverse_3d(out).unwrap();
        close(back.lam, 100.0, 1e-9, "tw3.inv.x");
        close(back.phi, 200.0, 1e-9, "tw3.inv.y");
        close(back.z, 300.0, 1e-9, "tw3.inv.z");
    }

    /// PROJ (verified against `cct +proj=helmert
    /// +towgs84=1,2,3,0.1,0.2,0.3,0.5`, no `+convention`) still requires an
    /// explicit `convention` even for `+towgs84` with non-zero rotations —
    /// error 1026 "missing 'convention' argument". `+towgs84` restricts
    /// *which* convention is accepted (see the next two tests); it does not
    /// exempt the caller from stating one.
    #[test]
    fn towgs84_seven_param_without_convention_is_missing_arg() {
        let ell = wgs84();
        let lk = LookupWithTowgs84::new("1.0,2.0,3.0,0.1,0.2,0.3,0.5");
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk,
            registry: None,
        };
        assert_eq!(new(&pp).err(), Some(ProjError::MissingArg));
    }

    /// PROJ (verified against `cct +proj=helmert
    /// +towgs84=1,2,3,0.1,0.2,0.3,0.5 +convention=coordinate_frame`) rejects
    /// `+towgs84` combined with an explicit `coordinate_frame` convention —
    /// error 1027 "towgs84 should only be used with
    /// convention=position_vector" — since the EPSG towgs84 shorthand is
    /// defined as position-vector (EPSG 9606) and PROJ hard-wires that.
    #[test]
    fn towgs84_seven_param_rejects_coordinate_frame() {
        let ell = wgs84();
        let lk = LookupWithTowgs84::new("1.0,2.0,3.0,0.1,0.2,0.3,0.5")
            .with_convention("coordinate_frame");
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk,
            registry: None,
        };
        assert_eq!(new(&pp).err(), Some(ProjError::IllegalArgValue));
    }

    /// PROJ (verified against `cct +proj=helmert
    /// +towgs84=1,2,3,0.1,0.2,0.3,0.5 +convention=position_vector`) accepts
    /// `+towgs84` combined with the matching `position_vector` convention.
    #[test]
    fn towgs84_seven_param_accepts_position_vector() {
        let ell = wgs84();
        let lk = LookupWithTowgs84::new("1.0,2.0,3.0,0.1,0.2,0.3,0.5")
            .with_convention("position_vector");
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk,
            registry: None,
        };
        assert!(new(&pp).is_ok());
    }

    #[test]
    fn transpose_flag_is_error() {
        let ell = wgs84();
        let lk = LookupWithTranspose::new();
        let pp = TransParams {
            ellipsoid: &ell,
            params: &lk,
            registry: None,
        };
        assert_eq!(new(&pp).err(), Some(ProjError::IllegalArgValue));
    }
}