oxiproj-transformations 0.1.2

//! Epoch-aware coordinate propagation using plate-motion models.
//!
//! Implements the foundation for T3.1 — automatic epoch-aware pathfinding.
//! Provides [`propagate_epoch`] for advancing ECEF coordinates through time
//! using a specified plate-motion model, and [`plate_for_frame`] for
//! heuristic frame-name to plate association.

use crate::frame_chain::{apply_frame_path, find_frame_path};
use oxiproj_core::{
    epoch::Epoch,
    plate_motion::{Plate, PlateMotionModel},
    ProjError, ProjResult,
};

/// Advance ECEF coordinates from one epoch to another using a plate-motion model.
///
/// Applies the velocity field encoded in `model` to the given geocentric
/// Cartesian coordinates `(x, y, z)`, propagating them from `from_epoch` to
/// `to_epoch`. The time delta is derived from the difference in decimal years
/// between the two epochs.
///
/// # Arguments
///
/// * `x` — ECEF X coordinate in metres.
/// * `y` — ECEF Y coordinate in metres.
/// * `z` — ECEF Z coordinate in metres.
/// * `from_epoch` — The epoch at which `(x, y, z)` is expressed.
/// * `to_epoch` — The target epoch to propagate the coordinate to.
/// * `model` — The plate-motion model to use for the propagation.
/// * `plate` — The tectonic plate on which the point lies.
///
/// # Errors
///
/// Returns [`ProjError::NoOperation`] when `model` contains no velocity
/// parameters for the specified `plate` (i.e., the model returns `None`).
///
/// # Examples
///
/// ```no_run
/// use oxiproj_core::epoch::Epoch;
/// use oxiproj_core::plate_motion::{Plate, PlateMotionModel};
/// use oxiproj_transformations::epoch_transform::propagate_epoch;
///
/// let from = Epoch::from_decimal_year(2015.0).unwrap();
/// let to   = Epoch::from_decimal_year(2025.0).unwrap();
/// let (x2, y2, z2) = propagate_epoch(
///     3_900_000.0, 900_000.0, 4_900_000.0,
///     &from, &to,
///     PlateMotionModel::Itrf2014, Plate::Eurasian,
/// ).unwrap();
/// ```
pub fn propagate_epoch(
    x: f64,
    y: f64,
    z: f64,
    from_epoch: &Epoch,
    to_epoch: &Epoch,
    model: PlateMotionModel,
    plate: Plate,
) -> ProjResult<(f64, f64, f64)> {
    match model.propagate(plate, x, y, z, from_epoch.year, to_epoch.year) {
        Some(result) => Ok(result),
        None => Err(ProjError::NoOperation),
    }
}

/// Heuristically map a reference-frame name string to a [`Plate`].
///
/// Performs a case-insensitive substring search on `frame_name` to identify
/// which tectonic plate a given reference-frame string is associated with.
/// This is a best-effort lookup intended for PROJ pipeline initialisation and
/// epoch-aware path selection; it should not be used as a substitute for a
/// proper authority-database lookup.
///
/// Returns `None` when no pattern matches.
///
/// # Pattern priority
///
/// Patterns are evaluated in the order listed below; the first match wins.
///
/// | Pattern | Returns |
/// |---------|---------|
/// | `"noam"` or `"north_am"` | [`Plate::NorthAmerican`] |
/// | `"euras"` | [`Plate::Eurasian`] |
/// | `"afr"` | [`Plate::African`] |
/// | `"pacific"` | [`Plate::Pacific`] |
/// | `"antarcti"` | [`Plate::Antarctic`] |
/// | `"austral"` | [`Plate::Australian`] |
/// | `"south_am"` | [`Plate::SouthAmerican`] |
/// | `"sam"` | [`Plate::SouthAmerican`] |
/// | `"arabian"` or `"arab"` | [`Plate::Arabian`] |
///
/// # Examples
///
/// ```
/// use oxiproj_core::plate_motion::Plate;
/// use oxiproj_transformations::epoch_transform::plate_for_frame;
///
/// assert_eq!(plate_for_frame("NOAM_2014"), Some(Plate::NorthAmerican));
/// assert_eq!(plate_for_frame("EURAS_plate"), Some(Plate::Eurasian));
/// assert_eq!(plate_for_frame("unknown_frame"), None);
/// ```
pub fn plate_for_frame(frame_name: &str) -> Option<Plate> {
    let lower = frame_name.to_ascii_lowercase();

    if lower.contains("noam") || lower.contains("north_am") {
        return Some(Plate::NorthAmerican);
    }
    if lower.contains("euras") {
        return Some(Plate::Eurasian);
    }
    if lower.contains("afr") {
        return Some(Plate::African);
    }
    if lower.contains("pacific") {
        return Some(Plate::Pacific);
    }
    if lower.contains("antarcti") {
        return Some(Plate::Antarctic);
    }
    if lower.contains("austral") {
        return Some(Plate::Australian);
    }
    // "south_am" must be tested before bare "sam" to avoid false positives
    if lower.contains("south_am") {
        return Some(Plate::SouthAmerican);
    }
    if lower.contains("sam") {
        return Some(Plate::SouthAmerican);
    }
    // "arabian" before "arab" — "arabian" already contains "arab", so either
    // order is equivalent here, but listing the longer string first is clearer.
    if lower.contains("arabian") || lower.contains("arab") {
        return Some(Plate::Arabian);
    }

    None
}

/// Multi-step, multi-hop epoch-aware pathfinding:
/// 1. Propagate from `from_epoch` to `to_epoch` in `from_frame` using plate motion
///    (the point physically moves on its plate; the frame is unchanged).
/// 2. Apply the frame-change Helmert **chain** if frames differ — a
///    breadth-first path through [`crate::frame_chain::FRAME_TRANSFORMS`],
///    evaluated at `to_epoch` (frame changes never alter the coordinate epoch).
/// 3. Result is in `to_frame` at `to_epoch`.
///
/// This resolves multi-hop cases (e.g. `ITRF2020→ITRF2000`, which composes two
/// hops through the ITRF2020 hub) that a single direct/inverse lookup cannot.
///
/// # Arguments
/// * `from_frame` — Source frame name (e.g. `"ITRF2020"`)
/// * `to_frame` — Target frame name (e.g. `"ITRF2000"`)
/// * `from_epoch` — Source epoch
/// * `to_epoch` — Target epoch
/// * `x, y, z` — ECEF coordinates in metres
/// * `model` — Plate motion model
/// * `plate` — Tectonic plate the point lies on (for the propagation step)
///
/// # Errors
///
/// Returns [`ProjError::NoOperation`] when the plate is absent from `model`, or
/// when no path of frame transforms connects `from_frame` to `to_frame`.
#[allow(clippy::too_many_arguments)]
pub fn epoch_path(
    from_frame: &str,
    to_frame: &str,
    from_epoch: &Epoch,
    to_epoch: &Epoch,
    x: f64,
    y: f64,
    z: f64,
    model: PlateMotionModel,
    plate: Plate,
) -> ProjResult<(f64, f64, f64)> {
    let (x1, y1, z1) = propagate_epoch(x, y, z, from_epoch, to_epoch, model, plate)?;

    let path = find_frame_path(from_frame, to_frame).ok_or(ProjError::NoOperation)?;
    apply_frame_path(&path, x1, y1, z1, to_epoch.year).map_err(|_| ProjError::NoOperation)
}

/// Ellipsoid-independent geocentric-to-geographic longitude/latitude of an ECEF
/// point, used purely for coarse plate identification.
///
/// Returns `(lon_deg, lat_deg)` where `lon` is the exact geodetic longitude
/// (`atan2(y, x)`) and `lat` is the **geocentric** latitude
/// (`atan2(z, hypot(x, y))`). Geocentric vs geodetic latitude differ by at most
/// ~0.19°, far below the resolution of the plate-boundary polygons, so no
/// ellipsoid parameters are required. Returns `None` if `(x, y, z)` is the
/// geocentre (degenerate).
fn ecef_to_lonlat_deg(x: f64, y: f64, z: f64) -> Option<(f64, f64)> {
    let p = (x * x + y * y).sqrt();
    if p == 0.0 && z == 0.0 {
        return None;
    }
    let lon = y.atan2(x).to_degrees();
    let lat = z.atan2(p).to_degrees();
    Some((lon, lat))
}

/// Coarse plate-boundary polygon (a closed ring of `(lon_deg, lat_deg)`
/// vertices) associated with a [`Plate`].
struct PlatePolygon {
    plate: Plate,
    ring: &'static [(f64, f64)],
}

/// Test whether `(lon, lat)` lies inside a geographic polygon ring.
///
/// Uses even-odd ray casting. Longitudes are unwrapped to within ±180° of the
/// ring's first vertex before the test, so rings that straddle the
/// antimeridian (e.g. the Pacific) are handled correctly for spans below 360°.
fn point_in_ring(lon: f64, lat: f64, ring: &[(f64, f64)]) -> bool {
    if ring.len() < 3 {
        return false;
    }
    let reference = ring[0].0;
    let unwrap = |l: f64| -> f64 {
        let mut v = l;
        while v - reference > 180.0 {
            v -= 360.0;
        }
        while v - reference < -180.0 {
            v += 360.0;
        }
        v
    };

    let plon = unwrap(lon);
    let mut inside = false;
    let n = ring.len();
    let mut j = n - 1;
    for i in 0..n {
        let (xi, yi) = (unwrap(ring[i].0), ring[i].1);
        let (xj, yj) = (unwrap(ring[j].0), ring[j].1);
        let intersects =
            ((yi > lat) != (yj > lat)) && (plon < (xj - xi) * (lat - yi) / (yj - yi) + xi);
        if intersects {
            inside = !inside;
        }
        j = i;
    }
    inside
}

/// Identify the tectonic [`Plate`] a geographic location lies on.
///
/// Performs a point-in-polygon lookup over a bundled set of coarse plate
/// boundaries covering the plates carried by the embedded plate-motion models.
/// Antarctica is handled by a latitude threshold (everything south of ~60°S).
/// Polygons are tested in priority order and the first containing ring wins.
///
/// Returns `None` for locations not covered by any polygon (chiefly open-ocean
/// interiors of plates that have no velocity in the models) — callers should
/// treat this as "plate unknown" rather than silently substituting a default.
///
/// The boundaries are deliberately coarse (continental-scale); they are meant
/// to select the correct rigid-plate velocity for epoch propagation, not to
/// resolve positions near a plate boundary to the metre.
#[must_use]
pub fn plate_at_lonlat(lon_deg: f64, lat_deg: f64) -> Option<Plate> {
    if lat_deg <= -60.0 {
        return Some(Plate::Antarctic);
    }
    for poly in PLATE_POLYGONS {
        if point_in_ring(lon_deg, lat_deg, poly.ring) {
            return Some(poly.plate);
        }
    }
    None
}

/// Identify the tectonic [`Plate`] an ECEF geocentric coordinate lies on.
///
/// Convenience wrapper over [`plate_at_lonlat`] that first converts the
/// geocentric `(x, y, z)` (metres) to a coarse longitude/latitude via
/// `ecef_to_lonlat_deg`. Returns `None` at the geocentre or for uncovered
/// locations.
#[must_use]
pub fn plate_for_coord(x: f64, y: f64, z: f64) -> Option<Plate> {
    let (lon, lat) = ecef_to_lonlat_deg(x, y, z)?;
    plate_at_lonlat(lon, lat)
}

/// Bundled coarse plate-boundary polygons, in lookup-priority order.
///
/// Rings are continental-scale approximations of the major tectonic plates
/// (longitude/latitude in degrees, closed rings). They are sufficient to select
/// the correct rigid-plate rotation for epoch propagation; they are **not**
/// authoritative plate boundaries. Antarctica is handled separately by a
/// latitude threshold in [`plate_at_lonlat`].
static PLATE_POLYGONS: &[PlatePolygon] = &[
    // Arabian plate — tested before Eurasian/African so it wins its region.
    PlatePolygon {
        plate: Plate::Arabian,
        ring: &[
            (34.0, 12.0),
            (34.0, 30.0),
            (42.0, 37.0),
            (48.0, 30.0),
            (58.0, 25.0),
            (60.0, 20.0),
            (52.0, 13.0),
            (43.0, 12.0),
            (40.0, 15.0),
        ],
    },
    // Eurasian plate (Europe + most of Asia north of the Alpine-Himalayan belt).
    PlatePolygon {
        plate: Plate::Eurasian,
        ring: &[
            (-11.0, 36.0),
            (-11.0, 44.0),
            (-10.0, 52.0),
            (2.0, 60.0),
            (5.0, 71.0),
            (30.0, 73.0),
            (60.0, 78.0),
            (100.0, 78.0),
            (140.0, 74.0),
            (145.0, 60.0),
            (145.0, 52.0),
            (135.0, 44.0),
            (122.0, 40.0),
            (108.0, 34.0),
            (98.0, 34.0),
            (88.0, 30.0),
            (75.0, 33.0),
            (62.0, 40.0),
            (44.0, 39.0),
            (40.0, 36.0),
            (28.0, 36.0),
            (10.0, 38.0),
        ],
    },
    // North American plate (North America + Greenland).
    PlatePolygon {
        plate: Plate::NorthAmerican,
        ring: &[
            (-168.0, 52.0),
            (-168.0, 65.0),
            (-150.0, 72.0),
            (-90.0, 82.0),
            (-20.0, 83.0),
            (-10.0, 62.0),
            (-55.0, 50.0),
            (-60.0, 47.0),
            (-82.0, 24.0),
            (-98.0, 18.0),
            (-106.0, 23.0),
            (-115.0, 30.0),
            (-125.0, 40.0),
            (-130.0, 50.0),
            (-140.0, 58.0),
        ],
    },
    // South American plate.
    PlatePolygon {
        plate: Plate::SouthAmerican,
        ring: &[
            (-82.0, 12.0),
            (-60.0, 12.0),
            (-50.0, 5.0),
            (-35.0, -5.0),
            (-34.0, -23.0),
            (-48.0, -34.0),
            (-58.0, -40.0),
            (-66.0, -55.0),
            (-75.0, -52.0),
            (-73.0, -30.0),
            (-70.0, -18.0),
            (-81.0, -6.0),
            (-82.0, 5.0),
        ],
    },
    // African plate (Nubia + Somalia, into the eastern Atlantic).
    PlatePolygon {
        plate: Plate::African,
        ring: &[
            (-18.0, 34.0),
            (10.0, 37.0),
            (24.0, 33.0),
            (34.0, 31.0),
            (43.0, 11.0),
            (51.0, 11.0),
            (42.0, -5.0),
            (40.0, -25.0),
            (35.0, -35.0),
            (18.0, -35.0),
            (12.0, -18.0),
            (8.0, 4.0),
            (-8.0, 5.0),
            (-17.0, 15.0),
        ],
    },
    // Australian plate (Australia + margin; Indo-Australian).
    PlatePolygon {
        plate: Plate::Australian,
        ring: &[
            (112.0, -10.0),
            (130.0, -9.0),
            (142.0, -9.0),
            (154.0, -24.0),
            (152.0, -40.0),
            (146.0, -45.0),
            (135.0, -45.0),
            (129.0, -38.0),
            (115.0, -37.0),
            (113.0, -22.0),
        ],
    },
    // Pacific plate — spans the antimeridian; ring longitudes run continuously
    // eastward from the western margin (Kamchatka/Japan trench) across 180° to
    // the East Pacific Rise, unwrapped relative to the first vertex.
    PlatePolygon {
        plate: Plate::Pacific,
        ring: &[
            (160.0, 55.0),
            (200.0, 52.0),
            (245.0, 40.0),
            (253.0, 24.0),
            (255.0, 5.0),
            (250.0, -20.0),
            (210.0, -55.0),
            (180.0, -30.0),
            (155.0, -5.0),
            (142.0, 20.0),
            (150.0, 40.0),
        ],
    },
];

#[cfg(test)]
mod tests {
    use super::*;
    use oxiproj_core::epoch::Epoch;
    use oxiproj_core::plate_motion::{Plate, PlateMotionModel};

    #[test]
    fn propagate_eurasian_plate_10_years() {
        let (x, y, z) = (3_900_000.0_f64, 900_000.0_f64, 4_900_000.0_f64);
        let from = Epoch::from_decimal_year(2015.0).unwrap();
        let to = Epoch::from_decimal_year(2025.0).unwrap();
        let (x2, y2, z2) = propagate_epoch(
            x,
            y,
            z,
            &from,
            &to,
            PlateMotionModel::Itrf2014,
            Plate::Eurasian,
        )
        .unwrap();
        let shift = ((x2 - x).powi(2) + (y2 - y).powi(2) + (z2 - z).powi(2)).sqrt();
        assert!(shift > 0.05 && shift < 0.5, "shift={}", shift);
    }

    #[test]
    fn plate_for_frame_lookups() {
        assert_eq!(plate_for_frame("NOAM_2014"), Some(Plate::NorthAmerican));
        assert_eq!(plate_for_frame("EURAS_plate"), Some(Plate::Eurasian));
        assert_eq!(plate_for_frame("unknown_frame"), None);
    }

    #[test]
    fn plate_for_frame_all_patterns() {
        assert_eq!(plate_for_frame("noam_itrf"), Some(Plate::NorthAmerican));
        assert_eq!(
            plate_for_frame("north_am_plate"),
            Some(Plate::NorthAmerican)
        );
        assert_eq!(plate_for_frame("AFR_2014"), Some(Plate::African));
        assert_eq!(plate_for_frame("PACIFIC"), Some(Plate::Pacific));
        assert_eq!(plate_for_frame("ANTARCTI_plate"), Some(Plate::Antarctic));
        assert_eq!(plate_for_frame("AUSTRAL_zone"), Some(Plate::Australian));
        assert_eq!(plate_for_frame("SOUTH_AM_2020"), Some(Plate::SouthAmerican));
        assert_eq!(plate_for_frame("ARABIAN_plate"), Some(Plate::Arabian));
        assert_eq!(plate_for_frame("ARAB_zone"), Some(Plate::Arabian));
        assert_eq!(plate_for_frame("UNKNOWN"), None);
    }

    #[test]
    fn propagate_returns_error_for_unknown_plate_in_model() {
        // NnrMorvel56 has no Arabian plate
        let from = Epoch::from_decimal_year(2015.0).unwrap();
        let to = Epoch::from_decimal_year(2025.0).unwrap();
        let result = propagate_epoch(
            1_000_000.0,
            2_000_000.0,
            3_000_000.0,
            &from,
            &to,
            PlateMotionModel::NnrMorvel56,
            Plate::Arabian,
        );
        assert!(result.is_err());
    }

    #[test]
    fn epoch_path_same_frame_is_propagation_only() {
        // Same frame, 1-year propagation: should match direct propagation
        let from = Epoch::from_decimal_year(2020.0).unwrap();
        let to = Epoch::from_decimal_year(2021.0).unwrap();
        let (x, y, z) = (4_627_798.0_f64, 119_795.0_f64, 4_369_668.0_f64);
        let result = epoch_path(
            "ITRF2020",
            "ITRF2020",
            &from,
            &to,
            x,
            y,
            z,
            PlateMotionModel::Itrf2020,
            Plate::Eurasian,
        );
        assert!(
            result.is_ok(),
            "epoch_path same frame failed: {:?}",
            result.err()
        );
        let (x2, y2, z2) = result.unwrap();
        // Should be close to propagate_epoch result
        let direct = propagate_epoch(
            x,
            y,
            z,
            &from,
            &to,
            PlateMotionModel::Itrf2020,
            Plate::Eurasian,
        )
        .unwrap();
        assert!((x2 - direct.0).abs() < 1e-6, "x mismatch");
        assert!((y2 - direct.1).abs() < 1e-6, "y mismatch");
        assert!((z2 - direct.2).abs() < 1e-6, "z mismatch");
    }

    #[test]
    fn epoch_path_itrf2020_to_itrf2014() {
        let from = Epoch::from_decimal_year(2020.0).unwrap();
        let to = Epoch::from_decimal_year(2020.0).unwrap();
        let (x, y, z) = (4_627_798.0_f64, 119_795.0_f64, 4_369_668.0_f64);
        let result = epoch_path(
            "ITRF2020",
            "ITRF2014",
            &from,
            &to,
            x,
            y,
            z,
            PlateMotionModel::Itrf2020,
            Plate::Eurasian,
        );
        assert!(
            result.is_ok(),
            "epoch_path cross frame failed: {:?}",
            result.err()
        );
        let (xo, yo, zo) = result.unwrap();
        // Frame transform should shift coordinates by mm-level
        let dist = ((xo - x).powi(2) + (yo - y).powi(2) + (zo - z).powi(2)).sqrt();
        assert!(dist < 1.0, "frame shift > 1m: {}", dist); // should be mm-level
        assert!(dist > 0.0, "no shift at all");
    }

    #[test]
    fn plate_at_lonlat_continents() {
        // Central Europe → Eurasian
        assert_eq!(plate_at_lonlat(13.4, 52.5), Some(Plate::Eurasian)); // Berlin
                                                                        // Denver → North American
        assert_eq!(plate_at_lonlat(-105.0, 39.7), Some(Plate::NorthAmerican));
        // Central Australia → Australian
        assert_eq!(plate_at_lonlat(133.0, -25.0), Some(Plate::Australian));
        // Central Brazil → South American
        assert_eq!(plate_at_lonlat(-55.0, -12.0), Some(Plate::SouthAmerican));
        // Sahara / central Africa → African
        assert_eq!(plate_at_lonlat(18.0, 12.0), Some(Plate::African));
        // Central Saudi Arabia → Arabian
        assert_eq!(plate_at_lonlat(45.0, 24.0), Some(Plate::Arabian));
        // Antarctica (latitude threshold) → Antarctic
        assert_eq!(plate_at_lonlat(0.0, -80.0), Some(Plate::Antarctic));
        // Mid Pacific (Hawaii-ish, ~-155°E = 205°) → Pacific
        assert_eq!(plate_at_lonlat(-155.0, 20.0), Some(Plate::Pacific));
    }

    #[test]
    fn plate_at_lonlat_open_ocean_is_none() {
        // Mid-Atlantic, far from any continental polygon → no plate assigned.
        assert_eq!(plate_at_lonlat(-30.0, 10.0), None);
    }

    #[test]
    fn plate_for_coord_from_ecef() {
        // Denver-ish ECEF (~-105°E, ~40°N) resolves to North American.
        assert_eq!(
            plate_for_coord(-1_274_000.0, -4_733_000.0, 4_074_000.0),
            Some(Plate::NorthAmerican)
        );
        // Central-European ECEF resolves to Eurasian.
        assert_eq!(
            plate_for_coord(3_711_000.0, 1_023_000.0, 5_036_000.0),
            Some(Plate::Eurasian)
        );
        // Geocentre is degenerate.
        assert_eq!(plate_for_coord(0.0, 0.0, 0.0), None);
    }

    #[test]
    fn epoch_path_multi_hop_itrf2020_to_itrf2000() {
        // ITRF2020→ITRF2000 needs 2 hops (via itself as hub is direct actually);
        // ITRF2014→ITRF2000 needs the ITRF2020 hub. Verify the 2-hop chain runs
        // and matches the PROJ pipeline value:
        //   cct +proj=pipeline +step +inv +init=ITRF2020:ITRF2014
        //                      +step +init=ITRF2020:ITRF2000
        //   in : 4627798.0 119795.0 4369668.0 @2015 → 4627798.013556 119795.002020 4369667.976067
        let e = Epoch::from_decimal_year(2015.0).unwrap();
        let (x, y, z) = (4_627_798.0_f64, 119_795.0_f64, 4_369_668.0_f64);
        let (xo, yo, zo) = epoch_path(
            "ITRF2014",
            "ITRF2000",
            &e,
            &e,
            x,
            y,
            z,
            PlateMotionModel::Itrf2020,
            Plate::Eurasian,
        )
        .expect("2-hop epoch_path should succeed");
        assert!((xo - 4_627_798.013_556).abs() < 1e-5, "x={xo}");
        assert!((yo - 119_795.002_020).abs() < 1e-5, "y={yo}");
        assert!((zo - 4_369_667.976_067).abs() < 1e-5, "z={zo}");
    }
}