WIP

Author: helgee

crates/lox-core/src/f64/consts.rs

@@ -54,3 +54,9 @@ pub const SECONDS_BETWEEN_MJD_AND_J2000: f64 = 4453444800.0;
 
 /// The number of seconds between the J1950 Epoch and the J2000 Epoch.
 pub const SECONDS_BETWEEN_J1950_AND_J2000: f64 = 1577880000.0;
+
+/*
+ * Physical constants
+ */
+
+pub const ROTATION_RATE_EARTH: f64 = 7.2921150e-5;

crates/lox-earth/src/frames.rs

@@ -0,0 +1,21 @@
+// SPDX-FileCopyrightText: 2025 Helge Eichhorn <[email protected]>
+//
+// SPDX-License-Identifier: MPL-2.0
+
+use lox_frames::transformations::Rotation;
+use lox_time::{Time, time_scales::Tt};
+
+use crate::precession::BiasPrecessionIau2000;
+
+pub mod iers1996;
+pub mod iers2003;
+pub mod iers2010;
+
+pub fn icrf_to_j2000(time: Time<Tt>) -> Rotation {
+    let bp = BiasPrecessionIau2000::new(time);
+    Rotation::new(bp.rb)
+}
+
+pub fn j2000_to_icrf(time: Time<Tt>) -> Rotation {
+    icrf_to_j2000(time).transpose()
+}

crates/lox-earth/src/frames/iers1996.rs

@@ -0,0 +1,57 @@
+// SPDX-FileCopyrightText: 2025 Helge Eichhorn <[email protected]>
+//
+// SPDX-License-Identifier: MPL-2.0
+
+use glam::DVec3;
+use lox_bodies::{Earth, RotationalElements};
+use lox_frames::transformations::Rotation;
+use lox_time::{
+    Time,
+    julian_dates::JulianDate,
+    time_scales::{Tdb, Tt, Ut1},
+};
+
+use crate::{
+    ecliptic::MeanObliquity, nutation::Nutation, precession::PrecessionIau1976,
+    rotation::GreenwichApparentSiderealTime,
+};
+
+pub fn j2000_to_mod(time: Time<Tdb>) -> Rotation {
+    Rotation::new(PrecessionIau1976::matrix(time))
+}
+
+pub fn mod_to_j2000(time: Time<Tdb>) -> Rotation {
+    j2000_to_mod(time).transpose()
+}
+
+pub fn mod_to_tod(time: Time<Tdb>, corrections: Option<Nutation>) -> Rotation {
+    let nut = Nutation::iau1980(time) + corrections.unwrap_or_default();
+    let epsa = MeanObliquity::iau1980(time.with_scale(Tt));
+    Rotation::new(nut.nutation_matrix(epsa))
+}
+
+pub fn tod_to_mod(time: Time<Tdb>, corrections: Option<Nutation>) -> Rotation {
+    mod_to_tod(time, corrections).transpose()
+}
+
+pub fn tod_to_pef(tdb: Time<Tdb>, ut1: Time<Ut1>) -> Rotation {
+    let gst = GreenwichApparentSiderealTime::iau1994(ut1);
+    let t = tdb.seconds_since_j2000();
+    Rotation::new(gst.0.rotation_z()).with_angular_velocity(DVec3::new(
+        0.0,
+        0.0,
+        Earth.rotation_rate(t),
+    ))
+}
+
+pub fn pef_to_tod(tdb: Time<Tdb>, ut1: Time<Ut1>) -> Rotation {
+    tod_to_pef(tdb, ut1).transpose()
+}
+
+pub fn pef_to_itrf() -> Rotation {
+    todo!()
+}
+
+pub fn itrf_to_pef() -> Rotation {
+    pef_to_itrf().transpose()
+}

crates/lox-earth/src/frames/iers2003.rs

@@ -0,0 +1,95 @@
+// SPDX-FileCopyrightText: 2025 Helge Eichhorn <[email protected]>
+//
+// SPDX-License-Identifier: MPL-2.0
+
+use glam::DVec3;
+use lox_core::{f64::consts::ROTATION_RATE_EARTH, units::Angle};
+use lox_frames::transformations::Rotation;
+use lox_time::{
+    Time,
+    time_scales::{Tdb, Tt, Ut1},
+};
+
+use crate::{
+    cip::CipCoords, ecliptic::MeanObliquity, nutation::Nutation, precession::BiasPrecessionIau2000,
+    rotation::GreenwichApparentSiderealTime,
+};
+
+pub fn icrf_to_mod(time: Time<Tt>) -> Rotation {
+    let bp = BiasPrecessionIau2000::new(time);
+    Rotation::new(bp.rpb)
+}
+
+pub fn mod_to_icrf(time: Time<Tt>) -> Rotation {
+    icrf_to_mod(time).transpose()
+}
+
+pub fn mod_to_tod_iau2000a(time: Time<Tt>, corrections: Option<CipCoords>) -> Rotation {
+    let CipCoords { x: dx, y: dy } = corrections.unwrap_or_default();
+    let nut = Nutation::iau2000a(time.with_scale(Tdb));
+    let epsa = MeanObliquity::iau1980(time);
+    let numat = nut.nutation_matrix(epsa);
+
+    let BiasPrecessionIau2000 { rpb, .. } = BiasPrecessionIau2000::new(time);
+    let rnpb = numat * rpb;
+    let v1 = DVec3::new(dx.as_f64(), dy.as_f64(), 0.0);
+    let v2 = rnpb * v1;
+    let corrections = Nutation {
+        longitude: Angle::new(v2.x / epsa.0.sin()),
+        obliquity: Angle::new(v2.y),
+    };
+
+    Rotation::new((nut + corrections).nutation_matrix(epsa))
+}
+
+pub fn tod_to_mod_iau2000a(time: Time<Tt>, corrections: Option<CipCoords>) -> Rotation {
+    mod_to_tod_iau2000a(time, corrections).transpose()
+}
+
+pub fn mod_to_tod_iau2000b(time: Time<Tt>, corrections: Option<CipCoords>) -> Rotation {
+    let CipCoords { x: dx, y: dy } = corrections.unwrap_or_default();
+    let nut = Nutation::iau2000b(time.with_scale(Tdb));
+    let epsa = MeanObliquity::iau1980(time);
+    let numat = nut.nutation_matrix(epsa);
+
+    let BiasPrecessionIau2000 { rpb, .. } = BiasPrecessionIau2000::new(time);
+    let rnpb = numat * rpb;
+    let v1 = DVec3::new(dx.as_f64(), dy.as_f64(), 0.0);
+    let v2 = rnpb * v1;
+    let corrections = Nutation {
+        longitude: Angle::new(v2.x / epsa.0.sin()),
+        obliquity: Angle::new(v2.y),
+    };
+
+    Rotation::new((nut + corrections).nutation_matrix(epsa))
+}
+
+pub fn tod_to_mod_iau2000b(time: Time<Tt>, corrections: Option<CipCoords>) -> Rotation {
+    mod_to_tod_iau2000b(time, corrections).transpose()
+}
+
+pub fn tod_to_pef_iau2000a(tt: Time<Tt>, ut1: Time<Ut1>) -> Rotation {
+    let gst = GreenwichApparentSiderealTime::iau2000a(tt, ut1);
+    Rotation::new(gst.0.rotation_z()).with_angular_velocity(DVec3::new(
+        0.0,
+        0.0,
+        ROTATION_RATE_EARTH,
+    ))
+}
+
+pub fn pef_to_tod_iau2000a(tt: Time<Tt>, ut1: Time<Ut1>) -> Rotation {
+    tod_to_pef_iau2000a(tt, ut1).transpose()
+}
+
+pub fn tod_to_pef_iau2000b(tt: Time<Tt>, ut1: Time<Ut1>) -> Rotation {
+    let gst = GreenwichApparentSiderealTime::iau2000b(tt, ut1);
+    Rotation::new(gst.0.rotation_z()).with_angular_velocity(DVec3::new(
+        0.0,
+        0.0,
+        ROTATION_RATE_EARTH,
+    ))
+}
+
+pub fn pef_to_tod_iau2000b(tt: Time<Tt>, ut1: Time<Ut1>) -> Rotation {
+    tod_to_pef_iau2000b(tt, ut1).transpose()
+}

crates/lox-earth/src/frames/iers2010.rs

@@ -0,0 +1,3 @@
+// SPDX-FileCopyrightText: 2025 Helge Eichhorn <[email protected]>
+//
+// SPDX-License-Identifier: MPL-2.0

crates/lox-earth/src/iers.rs

@@ -0,0 +1,87 @@
+// SPDX-FileCopyrightText: 2025 Helge Eichhorn <[email protected]>
+//
+// SPDX-License-Identifier: MPL-2.0
+
+mod sealed {
+    pub trait Sealed {}
+    impl Sealed for super::Iers1996 {}
+    impl Sealed for super::Iers2003 {}
+    impl Sealed for super::Iers2010 {}
+    impl Sealed for super::IersConvention {}
+}
+
+pub trait IersConventionId: sealed::Sealed {
+    fn id(&self) -> usize;
+}
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub struct Iers1996;
+
+impl IersConventionId for Iers1996 {
+    fn id(&self) -> usize {
+        0
+    }
+}
+
+#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
+pub enum Iau2000 {
+    #[default]
+    A = 1,
+    B = 2,
+}
+
+#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
+pub struct Iers2003(pub Iau2000);
+
+impl IersConventionId for Iers2003 {
+    fn id(&self) -> usize {
+        self.0 as usize
+    }
+}
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub struct Iers2010;
+
+impl IersConventionId for Iers2010 {
+    fn id(&self) -> usize {
+        3
+    }
+}
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum IersConvention {
+    Iers1996,
+    Iers2003(Iau2000),
+    Iers2010,
+}
+
+impl IersConventionId for IersConvention {
+    fn id(&self) -> usize {
+        match self {
+            IersConvention::Iers1996 => Iers1996.id(),
+            IersConvention::Iers2003(iau2000) => Iers2003(*iau2000).id(),
+            IersConvention::Iers2010 => Iers2010.id(),
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use rstest::rstest;
+
+    use super::*;
+
+    #[rstest]
+    #[case(Iers1996, 0)]
+    #[case(Iers2003(Iau2000::A), 1)]
+    #[case(Iers2003(Iau2000::B), 2)]
+    #[case(Iers2010, 3)]
+    #[case(IersConvention::Iers1996, 0)]
+    #[case(IersConvention::Iers2003(Iau2000::A), 1)]
+    #[case(IersConvention::Iers2003(Iau2000::B), 2)]
+    #[case(IersConvention::Iers2010, 3)]
+    fn test_iers_convention_id<T: IersConventionId>(#[case] iers: T, #[case] exp: usize) {
+        let act = iers.id();
+        assert_eq!(act, exp);
+    }
+}

crates/lox-earth/src/itrf.rs

@@ -13,7 +13,7 @@ use glam::{DMat3, DVec3};
 use lox_bodies::{Earth, RotationalElements};
 use lox_core::f64::consts::SECONDS_PER_DAY;
 use lox_frames::dynamic::DynTransformError;
-use lox_frames::transformations::TryTransform;
+use lox_frames::transformations::TryRotation;
 use lox_frames::transformations::rotations::Rotation;
 use lox_frames::{Cirf, DynFrame, Icrf, transform_provider};
 use lox_time::Time;
@@ -32,14 +32,14 @@ pub fn icrf_to_cirf(centuries: f64) -> Rotation {
     Rotation::new(m)
 }
 
-impl<T> TryTransform<Icrf, Cirf, T> for EopProvider
+impl<T> TryRotation<Icrf, Cirf, T> for EopProvider
 where
     T: TimeScale + Copy,
     Self: TryOffset<T, Tdb, Error = EopProviderError>,
 {
     type Error = EopProviderError;
 
-    fn try_transform(
+    fn try_rotation(
         &self,
         _origin: Icrf,
         _target: Cirf,
@@ -75,10 +75,10 @@ pub enum DynTransformEopError {
     Eop(#[from] EopProviderError),
 }
 
-impl TryTransform<DynFrame, DynFrame, DynTimeScale> for EopProvider {
+impl TryRotation<DynFrame, DynFrame, DynTimeScale> for EopProvider {
     type Error = DynTransformEopError;
 
-    fn try_transform(
+    fn try_rotation(
         &self,
         _origin: DynFrame,
         _target: DynFrame,

crates/lox-earth/src/lib.rs

@@ -10,7 +10,9 @@ pub mod coordinate_transformations;
 pub mod ecliptic;
 pub mod eop;
 pub mod ephemeris;
+pub mod frames;
 pub mod fundamental;
+pub mod iers;
 pub mod itrf;
 pub mod nutation;
 pub mod precession;

crates/lox-earth/src/nutation.rs

@@ -32,7 +32,7 @@ pub struct Nutation {
 }
 
 impl Nutation {
-    fn nutation_matrix(&self, mean_obliquity: MeanObliquity) -> DMat3 {
+    pub fn nutation_matrix(&self, mean_obliquity: MeanObliquity) -> DMat3 {
         let rot1 = mean_obliquity.0.rotation_x();
         let rot2 = (-self.longitude).rotation_z();
         let rot3 = (-(self.obliquity + mean_obliquity.0)).rotation_x();

crates/lox-earth/src/precession.rs

@@ -45,6 +45,11 @@ impl PrecessionIau1976 {
             * tas2r;
         Self { zeta, z, theta }
     }
+
+    pub fn matrix(time: Time<Tdb>) -> DMat3 {
+        let prec = Self::new(Time::default(), time);
+        (-prec.z).rotation_z() * prec.theta.rotation_y() * (-prec.zeta).rotation_z()
+    }
 }
 
 pub struct BiasPrecessionIau2000 {
@@ -53,7 +58,7 @@ pub struct BiasPrecessionIau2000 {
     /// Precession matrix.
     pub rp: DMat3,
     /// Bias-precession matrix.
-    pub rbp: DMat3,
+    pub rpb: DMat3,
 }
 
 impl BiasPrecessionIau2000 {
@@ -81,7 +86,7 @@ impl BiasPrecessionIau2000 {
 
         // Bias-precession matrix: GCRS to mean of date
         let rbp = rp * rb;
-        Self { rb, rp, rbp }
+        Self { rb, rp, rpb: rbp }
     }
 }
 
@@ -225,7 +230,7 @@ mod tests {
             0.999_999_928_598_468_3,
         ])
         .transpose();
-        let rbp_exp = DMat3::from_cols_array(&[
+        let rpb_exp = DMat3::from_cols_array(&[
             0.999_999_550_517_508_8,
             8.695_405_883_617_885e-4,
             3.779_734_722_239_007e-4,
@@ -239,7 +244,7 @@ mod tests {
         .transpose();
         assert_approx_eq!(bp.rb, rb_exp, atol <= 1e-12);
         assert_approx_eq!(bp.rp, rp_exp, atol <= 1e-12);
-        assert_approx_eq!(bp.rbp, rbp_exp, atol <= 1e-12);
+        assert_approx_eq!(bp.rpb, rpb_exp, atol <= 1e-12);
     }
 
     #[test]