oxiproj-core 0.1.2

//! Plate-motion parameter catalog for OxiProj.
//!
//! This module provides angular velocity parameters for major tectonic plates
//! under three published plate-motion models:
//!
//! - **ITRF2014 PMM** — Altamimi et al., 2017, *Journal of Geophysical Research*
//! - **ITRF2020 PMM** — Altamimi et al., 2023
//! - **NNR-MORVEL56** — DeMets et al., 2010, *Geophysical Journal International*
//!
//! Angular velocities are stored in milli-arcseconds per year (mas/yr) and
//! converted to radians per year internally for velocity computations.

use core::f64::consts::PI;

/// Conversion factor: 1 milli-arcsecond/yr → radians/yr.
///
/// Derived as `PI / (180.0 * 3_600_000.0)`.
const MAS_TO_RAD: f64 = PI / (180.0 * 3_600_000.0);

/// Angular velocity of a tectonic plate expressed in the ITRF Cartesian frame.
///
/// All three components are in milli-arcseconds per year (mas/yr).
/// The vector `(omega_x, omega_y, omega_z)` is the Euler pole describing
/// rigid rotation of the plate about the geocentre.
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct PlateAngularVelocity {
    /// X-component of the angular velocity in milli-arcseconds per year.
    pub omega_x: f64,
    /// Y-component of the angular velocity in milli-arcseconds per year.
    pub omega_y: f64,
    /// Z-component of the angular velocity in milli-arcseconds per year.
    pub omega_z: f64,
}

/// Major tectonic plates supported by the embedded plate-motion models.
///
/// Not every model contains parameters for every plate; see
/// [`PlateMotionModel::angular_velocity`] for details.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[non_exhaustive]
pub enum Plate {
    /// Eurasian plate.
    Eurasian,
    /// North American plate.
    NorthAmerican,
    /// Pacific plate.
    Pacific,
    /// African plate.
    African,
    /// Antarctic plate.
    Antarctic,
    /// Australian plate.
    Australian,
    /// South American plate.
    SouthAmerican,
    /// Arabian plate.
    Arabian,
    /// Caribbean plate.
    Caribbean,
    /// Cocos plate.
    Cocos,
    /// Nazca plate.
    Nazca,
    /// Philippine plate.
    Philippine,
}

/// Published plate-motion model used for angular velocity lookups.
///
/// Each variant corresponds to a peer-reviewed set of Euler-pole parameters.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[non_exhaustive]
pub enum PlateMotionModel {
    /// ITRF2014 Plate Motion Model (Altamimi et al., 2017, *JGR*).
    Itrf2014,
    /// ITRF2020 Plate Motion Model (Altamimi et al., 2023).
    Itrf2020,
    /// NNR-MORVEL56 (DeMets et al., 2010, *Geophysical Journal International*).
    NnrMorvel56,
}

impl PlateMotionModel {
    /// Return the angular velocity of `plate` in this model, or `None` if the
    /// plate is not included in the model's parameter set.
    ///
    /// Angular velocity components are in milli-arcseconds per year (mas/yr).
    #[must_use]
    pub fn angular_velocity(&self, plate: Plate) -> Option<PlateAngularVelocity> {
        match self {
            PlateMotionModel::Itrf2014 => itrf2014_angular_velocity(plate),
            PlateMotionModel::Itrf2020 => itrf2020_angular_velocity(plate),
            PlateMotionModel::NnrMorvel56 => nnr_morvel56_angular_velocity(plate),
        }
    }

    /// Compute the ITRF Cartesian velocity (mm/yr) of a point on `plate`.
    ///
    /// `x`, `y`, `z` are geocentric coordinates in **metres**.
    ///
    /// Returns `Some((vx, vy, vz))` in **mm/yr**, or `None` if the plate is
    /// not present in this model.
    ///
    /// The velocity is computed as the cross product `ω × r`, where `ω` is the
    /// angular velocity vector converted to rad/yr and `r = (x, y, z)`.
    #[must_use]
    pub fn velocity_at(&self, plate: Plate, x: f64, y: f64, z: f64) -> Option<(f64, f64, f64)> {
        let av = self.angular_velocity(plate)?;

        let wx = av.omega_x * MAS_TO_RAD;
        let wy = av.omega_y * MAS_TO_RAD;
        let wz = av.omega_z * MAS_TO_RAD;

        // Cross product v = ω × r (result in m/yr when r is in metres)
        let vx = wy * z - wz * y;
        let vy = wz * x - wx * z;
        let vz = wx * y - wy * x;

        // Convert m/yr → mm/yr
        Some((vx * 1000.0, vy * 1000.0, vz * 1000.0))
    }

    /// Compute the coordinate shift (mm) accumulated over `dt_years` years
    /// for a point on `plate`.
    ///
    /// `x`, `y`, `z` are geocentric coordinates in **metres**.
    ///
    /// Returns `Some((dx, dy, dz))` in **mm**, or `None` if the plate is not
    /// present in this model.
    #[must_use]
    pub fn coordinate_shift(
        &self,
        plate: Plate,
        x: f64,
        y: f64,
        z: f64,
        dt_years: f64,
    ) -> Option<(f64, f64, f64)> {
        let (vx, vy, vz) = self.velocity_at(plate, x, y, z)?;
        Some((vx * dt_years, vy * dt_years, vz * dt_years))
    }

    /// Propagate a geocentric coordinate from `from_epoch` to `to_epoch` using
    /// the rigid-plate motion model.
    ///
    /// `x`, `y`, `z` are geocentric coordinates in **metres** at `from_epoch`.
    ///
    /// Returns `Some((x', y', z'))` in **metres** at `to_epoch`, or `None` if
    /// the plate is not present in this model.
    ///
    /// The propagation is a simple linear (first-order) model: the plate is
    /// assumed to rotate at constant angular velocity.
    #[must_use]
    pub fn propagate(
        &self,
        plate: Plate,
        x: f64,
        y: f64,
        z: f64,
        from_epoch: f64,
        to_epoch: f64,
    ) -> Option<(f64, f64, f64)> {
        let dt = to_epoch - from_epoch;
        let (dx_mm, dy_mm, dz_mm) = self.coordinate_shift(plate, x, y, z, dt)?;

        // Convert mm → m
        let dx_m = dx_mm / 1000.0;
        let dy_m = dy_mm / 1000.0;
        let dz_m = dz_mm / 1000.0;

        Some((x + dx_m, y + dy_m, z + dz_m))
    }
}

// ---------------------------------------------------------------------------
// Private lookup functions — one per model
// ---------------------------------------------------------------------------

/// ITRF2014 Plate Motion Model parameters (Altamimi et al., 2017, JGR).
fn itrf2014_angular_velocity(plate: Plate) -> Option<PlateAngularVelocity> {
    let (ox, oy, oz) = match plate {
        Plate::Eurasian => (-0.085, -0.531, 0.770),
        Plate::NorthAmerican => (0.024, -0.694, -0.033),
        Plate::Pacific => (-0.409, 1.047, -2.169),
        Plate::African => (0.000, 0.310, 0.550),
        Plate::Antarctic => (-0.252, -0.302, 0.643),
        Plate::Australian => (1.504, 1.172, 1.228),
        Plate::SouthAmerican => (-0.270, -0.301, -0.140),
        Plate::Arabian => (1.154, -0.136, 1.444),
        Plate::Caribbean | Plate::Cocos | Plate::Nazca | Plate::Philippine => return None,
    };
    Some(PlateAngularVelocity {
        omega_x: ox,
        omega_y: oy,
        omega_z: oz,
    })
}

/// ITRF2020 Plate Motion Model parameters (Altamimi et al., 2023).
fn itrf2020_angular_velocity(plate: Plate) -> Option<PlateAngularVelocity> {
    let (ox, oy, oz) = match plate {
        Plate::Eurasian => (-0.083, -0.534, 0.750),
        Plate::NorthAmerican => (0.029, -0.698, -0.030),
        Plate::Pacific => (-0.411, 1.049, -2.162),
        Plate::African => (0.002, 0.316, 0.546),
        Plate::Antarctic => (-0.247, -0.294, 0.635),
        Plate::Australian => (1.511, 1.166, 1.218),
        Plate::SouthAmerican => (-0.265, -0.308, -0.132),
        Plate::Arabian => (1.156, -0.136, 1.451),
        Plate::Caribbean | Plate::Cocos | Plate::Nazca | Plate::Philippine => return None,
    };
    Some(PlateAngularVelocity {
        omega_x: ox,
        omega_y: oy,
        omega_z: oz,
    })
}

/// NNR-MORVEL56 parameters (DeMets et al., 2010, Geophysical Journal International).
fn nnr_morvel56_angular_velocity(plate: Plate) -> Option<PlateAngularVelocity> {
    let (ox, oy, oz) = match plate {
        Plate::Eurasian => (-0.083, -0.534, 0.750),
        Plate::NorthAmerican => (0.035, -0.662, -0.100),
        Plate::Pacific => (-0.444, 1.049, -2.093),
        Plate::African => (0.048, 0.280, 0.564),
        Plate::Antarctic => (-0.223, -0.310, 0.594),
        Plate::Australian => (1.504, 1.172, 1.228),
        Plate::SouthAmerican => (-0.270, -0.301, -0.140),
        Plate::Arabian | Plate::Caribbean | Plate::Cocos | Plate::Nazca | Plate::Philippine => {
            return None
        }
    };
    Some(PlateAngularVelocity {
        omega_x: ox,
        omega_y: oy,
        omega_z: oz,
    })
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

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

    #[test]
    fn pacific_plate_velocity_non_zero() {
        let omega = PlateMotionModel::Itrf2014
            .angular_velocity(Plate::Pacific)
            .unwrap();
        assert!(omega.omega_z.abs() > 0.0);
    }

    #[test]
    fn eurasian_velocity_at_origin_is_zero() {
        let (vx, vy, vz) = PlateMotionModel::Itrf2014
            .velocity_at(Plate::Eurasian, 0.0, 0.0, 0.0)
            .unwrap();
        assert_eq!(vx, 0.0);
        assert_eq!(vy, 0.0);
        assert_eq!(vz, 0.0);
    }

    #[test]
    fn north_american_propagation_makes_sense() {
        let (x, y, z) = (-2_000_000.0_f64, -4_500_000.0_f64, 4_000_000.0_f64);
        let result = PlateMotionModel::Itrf2014
            .propagate(Plate::NorthAmerican, x, y, z, 2010.0, 2020.0)
            .unwrap();
        let shift =
            ((result.0 - x).powi(2) + (result.1 - y).powi(2) + (result.2 - z).powi(2)).sqrt();
        assert!(
            shift < 0.5,
            "10-year plate motion shift > 0.5m is unrealistic: {}",
            shift
        );
        assert!(
            shift > 0.001,
            "10-year plate motion shift should be detectable"
        );
    }

    #[test]
    fn all_itrf2014_plates_have_velocity() {
        let plates = [
            Plate::Eurasian,
            Plate::NorthAmerican,
            Plate::Pacific,
            Plate::African,
            Plate::Antarctic,
            Plate::Australian,
            Plate::SouthAmerican,
            Plate::Arabian,
        ];
        for plate in plates {
            assert!(
                PlateMotionModel::Itrf2014.angular_velocity(plate).is_some(),
                "{:?} missing from ITRF2014",
                plate
            );
        }
    }

    #[test]
    fn itrf2020_eurasian_slightly_different_from_itrf2014() {
        let v14 = PlateMotionModel::Itrf2014
            .angular_velocity(Plate::Eurasian)
            .unwrap();
        let v20 = PlateMotionModel::Itrf2020
            .angular_velocity(Plate::Eurasian)
            .unwrap();
        // Different models, slightly different values
        assert_ne!(v14.omega_z, v20.omega_z);
    }
}