oxiproj-transformations 0.1.1

//! ITRF frame-chain pathfinding and epoch-aware Helmert application.
//!
//! Provides a static catalogue of 7-parameter Helmert transforms with rates
//! between successive ITRF realisations, and utilities for finding and
//! applying them (forward and inverse) at an arbitrary decimal-year epoch.

use oxiproj_core::ProjResult;

/// A 7-parameter Helmert transform with temporal rates between two ITRF frames.
///
/// All parameters follow the IERS conventions:
/// - translations in mm (rate in mm/yr)
/// - scale in parts per billion (rate in ppb/yr)
/// - rotations in milli-arc-seconds (rate in mas/yr)
#[derive(Debug, Clone, Copy)]
pub struct FrameTransform {
    /// Source frame name (e.g. `"ITRF2020"`).
    pub from: &'static str,
    /// Target frame name (e.g. `"ITRF2014"`).
    pub to: &'static str,
    /// Reference epoch (decimal year) at which the published parameters apply.
    pub ref_epoch: f64,
    // Translations (mm)
    pub tx_mm: f64,
    pub ty_mm: f64,
    pub tz_mm: f64,
    // Translation rates (mm/yr)
    pub dtx_mm: f64,
    pub dty_mm: f64,
    pub dtz_mm: f64,
    // Scale (ppb) and rate (ppb/yr)
    pub scale_ppb: f64,
    pub dscale_ppb: f64,
    // Rotations (mas) and rates (mas/yr)
    pub rx_mas: f64,
    pub ry_mas: f64,
    pub rz_mas: f64,
    pub drx_mas: f64,
    pub dry_mas: f64,
    pub drz_mas: f64,
}

/// Static catalogue of inter-ITRF Helmert transforms.
///
/// Parameters are taken from the IERS published transformation tables.
/// Each entry represents the transform FROM `from` TO `to` at `ref_epoch`.
pub const FRAME_TRANSFORMS: &[FrameTransform] = &[
    // ITRF2020 → ITRF2014  (ref_epoch = 2015.0)
    // Source: Altamimi et al. (2023), Table 2
    FrameTransform {
        from: "ITRF2020",
        to: "ITRF2014",
        ref_epoch: 2015.0,
        tx_mm: -1.4,
        ty_mm: 0.9,
        tz_mm: 1.4,
        dtx_mm: 0.0,
        dty_mm: -0.1,
        dtz_mm: 0.2,
        scale_ppb: -0.42,
        dscale_ppb: 0.0,
        rx_mas: 0.0,
        ry_mas: 0.0,
        rz_mas: 0.0,
        drx_mas: 0.0,
        dry_mas: 0.0,
        drz_mas: 0.0,
    },
    // ITRF2014 → ITRF2008  (ref_epoch = 2010.0)
    // Source: Altamimi et al. (2016), Table 2
    FrameTransform {
        from: "ITRF2014",
        to: "ITRF2008",
        ref_epoch: 2010.0,
        tx_mm: 1.6,
        ty_mm: 1.9,
        tz_mm: 2.4,
        dtx_mm: 0.0,
        dty_mm: 0.0,
        dtz_mm: -0.1,
        scale_ppb: -0.02,
        dscale_ppb: 0.03,
        rx_mas: 0.0,
        ry_mas: 0.0,
        rz_mas: 0.0,
        drx_mas: 0.0,
        dry_mas: 0.0,
        drz_mas: 0.0,
    },
];

/// Search [`FRAME_TRANSFORMS`] for a direct transform between two frames.
///
/// Comparison is case-insensitive. Returns `None` when no entry matches.
pub fn find_frame_transform<'a>(from: &str, to: &str) -> Option<&'a FrameTransform> {
    FRAME_TRANSFORMS
        .iter()
        .find(|ft| ft.from.eq_ignore_ascii_case(from) && ft.to.eq_ignore_ascii_case(to))
}

/// Apply a Helmert frame transform at a given decimal-year `epoch`.
///
/// Uses a linearised (small-angle) Helmert model:
///
/// ```text
/// [xout]   [tx]         [  1  -rz  ry ] [x]
/// [yout] = [ty] + scale·[ rz   1  -rx ]·[y]
/// [zout]   [tz]         [-ry  rx   1  ] [z]
/// ```
///
/// Parameters are propagated from `transform.ref_epoch` to `epoch` using the
/// published rates before the transform is applied.
///
/// # Errors
///
/// This function always succeeds; the return type is `ProjResult` for API
/// consistency with the rest of the pipeline.
pub fn apply_frame_transform(
    transform: &FrameTransform,
    x: f64,
    y: f64,
    z: f64,
    epoch: f64,
) -> ProjResult<(f64, f64, f64)> {
    let dt = epoch - transform.ref_epoch;

    let tx = (transform.tx_mm + transform.dtx_mm * dt) * 1e-3;
    let ty = (transform.ty_mm + transform.dty_mm * dt) * 1e-3;
    let tz = (transform.tz_mm + transform.dtz_mm * dt) * 1e-3;

    let scale = 1.0 + (transform.scale_ppb + transform.dscale_ppb * dt) * 1e-9;

    // 1 mas = π / (180 × 3_600_000) radians
    const MAS_TO_RAD: f64 = core::f64::consts::PI / (180.0 * 3_600_000.0);

    let rx = (transform.rx_mas + transform.drx_mas * dt) * MAS_TO_RAD;
    let ry = (transform.ry_mas + transform.dry_mas * dt) * MAS_TO_RAD;
    let rz = (transform.rz_mas + transform.drz_mas * dt) * MAS_TO_RAD;

    let xout = tx + scale * (x - rz * y + ry * z);
    let yout = ty + scale * (rz * x + y - rx * z);
    let zout = tz + scale * (-ry * x + rx * y + z);

    Ok((xout, yout, zout))
}

/// Apply the inverse of a Helmert frame transform at a given decimal-year `epoch`.
///
/// Constructs an inverted [`FrameTransform`] by negating all translation,
/// scale, and rotation parameters (and their rates), then delegates to
/// [`apply_frame_transform`].
pub fn apply_frame_transform_inverse(
    transform: &FrameTransform,
    x: f64,
    y: f64,
    z: f64,
    epoch: f64,
) -> ProjResult<(f64, f64, f64)> {
    let inv = FrameTransform {
        from: transform.from,
        to: transform.to,
        ref_epoch: transform.ref_epoch,
        tx_mm: -transform.tx_mm,
        ty_mm: -transform.ty_mm,
        tz_mm: -transform.tz_mm,
        dtx_mm: -transform.dtx_mm,
        dty_mm: -transform.dty_mm,
        dtz_mm: -transform.dtz_mm,
        scale_ppb: -transform.scale_ppb,
        dscale_ppb: -transform.dscale_ppb,
        rx_mas: -transform.rx_mas,
        ry_mas: -transform.ry_mas,
        rz_mas: -transform.rz_mas,
        drx_mas: -transform.drx_mas,
        dry_mas: -transform.dry_mas,
        drz_mas: -transform.drz_mas,
    };
    apply_frame_transform(&inv, x, y, z, epoch)
}

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

    #[test]
    fn frame_transform_itrf2020_to_itrf2014_found() {
        let ft = find_frame_transform("ITRF2020", "ITRF2014");
        assert!(ft.is_some(), "ITRF2020->ITRF2014 transform not found");
    }

    #[test]
    fn frame_transform_itrf2020_to_itrf2014_apply() {
        let ft = find_frame_transform("ITRF2020", "ITRF2014").unwrap();
        let result = apply_frame_transform(ft, 4_627_798.0, 119_795.0, 4_369_668.0, 2015.0);
        assert!(result.is_ok());
        let (x, y, z) = result.unwrap();
        // < 1m difference from input (mm-level transform)
        assert!((x - 4_627_798.0).abs() < 1.0);
        assert!((y - 119_795.0).abs() < 1.0);
        assert!((z - 4_369_668.0).abs() < 1.0);
    }

    #[test]
    fn frame_transform_inverse_reverses_forward() {
        let ft = find_frame_transform("ITRF2020", "ITRF2014").unwrap();
        let x0 = 4_627_798.0_f64;
        let y0 = 119_795.0_f64;
        let z0 = 4_369_668.0_f64;
        let (xf, yf, zf) = apply_frame_transform(ft, x0, y0, z0, 2015.0).unwrap();
        let (xr, yr, zr) = apply_frame_transform_inverse(ft, xf, yf, zf, 2015.0).unwrap();
        // Roundtrip should be sub-mm accurate
        assert!(
            (xr - x0).abs() < 0.001,
            "x roundtrip error: {}",
            (xr - x0).abs()
        );
        assert!(
            (yr - y0).abs() < 0.001,
            "y roundtrip error: {}",
            (yr - y0).abs()
        );
        assert!(
            (zr - z0).abs() < 0.001,
            "z roundtrip error: {}",
            (zr - z0).abs()
        );
    }

    #[test]
    fn frame_transform_unknown_returns_none() {
        assert!(find_frame_transform("ITRF1900", "ITRF2014").is_none());
    }
}