feat(lox-earth): implement celestial to terrestrial rotations

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/cio.rs

@@ -5,4 +5,10 @@
 
 //! Module cio exposes functions for calculating the Celestial Intermediate Origin (CIO) locator, s.
 
-pub mod s06;
+use lox_core::units::Angle;
+use lox_test_utils::ApproxEq;
+
+pub mod iau2006;
+
+#[derive(Debug, Clone, Copy, Default, PartialEq, ApproxEq)]
+pub struct CioLocator(pub Angle);

crates/lox-earth/src/cio/s06.rs renamed to crates/lox-earth/src/cio/iau2006.rs

@@ -6,9 +6,14 @@
 //! Module s06 exposes a function for calculating the Celestial Intermediate Origin (CIO) locator,
 //! s, using IAU 2006 precession and IAU 2000A nutation.
 
+use fast_polynomial::poly_array;
 use lox_core::types::units::JulianCenturies;
 use lox_core::units::{Angle, AngleUnits};
+use lox_time::Time;
+use lox_time::julian_dates::JulianDate;
+use lox_time::time_scales::Tdb;
 
+use crate::cio::CioLocator;
 use crate::cip::CipCoords;
 use crate::fundamental::iers03::{
     d_iers03, earth_l_iers03, f_iers03, l_iers03, lp_iers03, omega_iers03, pa_iers03,
@@ -17,6 +22,16 @@ use crate::fundamental::iers03::{
 
 mod terms;
 
+impl CioLocator {
+    pub fn iau2006(time: Time<Tdb>, CipCoords { x, y }: CipCoords) -> Self {
+        let t = time.centuries_since_j2000();
+        let fundamental_args = fundamental_args(t);
+        let evaluated_terms = evaluate_terms(&fundamental_args);
+        let s = poly_array(t, &evaluated_terms).arcsec();
+        CioLocator(s - y.to_radians() / 2.0 * x)
+    }
+}
+
 /// l, l', F, D, Ω, LVe, LE and pA.
 type FundamentalArgs = [Angle; 8];
 
@@ -28,7 +43,7 @@ pub fn cio_locator(
 ) -> Angle {
     let fundamental_args = fundamental_args(centuries_since_j2000_tdb);
     let evaluated_terms = evaluate_terms(&fundamental_args);
-    let s = fast_polynomial::poly_array(centuries_since_j2000_tdb, &evaluated_terms).arcsec();
+    let s = poly_array(centuries_since_j2000_tdb, &evaluated_terms).arcsec();
     Angle::radians(s.to_radians() - x.to_radians() * y.to_radians() / 2.0)
 }
 
@@ -76,52 +91,27 @@ fn evaluate_single_order_terms(
 
 #[cfg(test)]
 mod tests {
-    use std::iter::zip;
-
     use lox_test_utils::assert_approx_eq;
 
     use super::*;
 
-    const TOLERANCE: f64 = 1e-11;
-
-    #[test]
-    fn test_s_jd0() {
-        let jd0: JulianCenturies = -67.11964407939767;
-        let xy = CipCoords::new(jd0);
-        assert_approx_eq!(
-            cio_locator(jd0, xy),
-            -0.0723985415686306.rad(),
-            rtol <= TOLERANCE
-        );
-    }
-
-    #[test]
-    fn test_s_j2000() {
-        let j2000: JulianCenturies = 0.0;
-        let xy = CipCoords::new(j2000);
-        assert_approx_eq!(
-            cio_locator(j2000, xy),
-            -0.00000001013396519178.rad(),
-            rtol <= TOLERANCE
-        );
-    }
-
     #[test]
-    fn test_s_j2100() {
-        let j2100: JulianCenturies = 1.0;
-        let xy = CipCoords::new(j2100);
-        assert_approx_eq!(
-            cio_locator(j2100, xy),
-            -0.00000000480511934533.rad(),
-            rtol <= TOLERANCE
-        );
+    fn test_cio_locator_iau2006() {
+        let cip_coords = CipCoords {
+            x: 5.791_308_486_706_011e-4.rad(),
+            y: 4.020_579_816_732_961e-5.rad(),
+        };
+        let time = Time::from_two_part_julian_date(Tdb, 2400000.5, 53736.0);
+        let exp = CioLocator(-1.220_032_213_076_463e-8.rad());
+        let act = CioLocator::iau2006(time, cip_coords);
+        assert_approx_eq!(act, exp, atol <= 1e-18);
     }
 
     #[test]
     fn test_fundamental_args_ordering() {
         let j2000: JulianCenturies = 0.0;
-        let actual = fundamental_args(j2000);
-        let expected = [
+        let act = fundamental_args(j2000);
+        let exp = [
             l_iers03(j2000),
             lp_iers03(j2000),
             f_iers03(j2000),
@@ -132,8 +122,6 @@ mod tests {
             pa_iers03(j2000),
         ];
 
-        for (act, exp) in zip(actual, expected) {
-            assert_approx_eq!(act, exp, rtol <= TOLERANCE)
-        }
+        assert_approx_eq!(act, exp, rtol <= 1e-12)
     }
 }

crates/lox-earth/src/cio/s06/terms.rs renamed to crates/lox-earth/src/cio/iau2006/terms.rs

@@ -6,13 +6,15 @@
 //! Module cip exposes functions for calculating the position of the
 //! Celestial Intermediate Pole (CIP).
 
+use std::ops::{Add, AddAssign};
+
 use glam::DMat3;
-use lox_core::{
-    types::units::JulianCenturies,
-    units::{Angle, AngleUnits},
-};
+use lox_core::units::{Angle, AngleUnits};
+use lox_time::{Time, julian_dates::JulianDate, time_scales::Tdb};
+
+use crate::cio::CioLocator;
 
-mod xy06;
+mod iau2006;
 
 #[derive(Copy, Clone, Debug, Default, PartialEq)]
 pub struct CipCoords {
@@ -23,8 +25,9 @@ pub struct CipCoords {
 impl CipCoords {
     /// Calculates the (X, Y) coordinates of the Celestial Intermediate Pole (CIP) using the the IAU
     /// 2006 precession and IAU 2000A nutation models.
-    pub fn new(centuries_since_j2000_tdb: JulianCenturies) -> Self {
-        let (x, y) = xy06::cip_coords(centuries_since_j2000_tdb);
+    pub fn iau2006(time: Time<Tdb>) -> Self {
+        let t = time.centuries_since_j2000();
+        let (x, y) = iau2006::cip_coords(t);
         Self { x, y }
     }
 
@@ -34,10 +37,45 @@ impl CipCoords {
         let y = bpn.y_axis.z.rad();
         Self { x, y }
     }
+
+    pub fn celestial_to_intermediate_matrix(&self, s: CioLocator) -> DMat3 {
+        let x = self.x.to_radians();
+        let y = self.y.to_radians();
+        let r2 = x.powi(2) + y.powi(2);
+
+        let e = if r2 > 0.0 {
+            Angle::from_atan2(y, x)
+        } else {
+            Angle::default()
+        };
+        let d = Angle::from_atan((r2 / (1.0 - r2)).sqrt());
+
+        (-(e + s.0)).rotation_z() * d.rotation_y() * e.rotation_z()
+    }
+}
+
+impl Add for CipCoords {
+    type Output = Self;
+
+    fn add(self, rhs: Self) -> Self::Output {
+        CipCoords {
+            x: self.x + rhs.x,
+            y: self.y + rhs.y,
+        }
+    }
+}
+
+impl AddAssign for CipCoords {
+    fn add_assign(&mut self, rhs: Self) {
+        self.x += rhs.x;
+        self.y += rhs.y;
+    }
 }
 
 #[cfg(test)]
 mod tests {
+    use lox_test_utils::assert_approx_eq;
+
     use super::*;
 
     #[test]
@@ -58,4 +96,27 @@ mod tests {
         assert_eq!(cip.x, 1.093465510215479e-3.rad());
         assert_eq!(cip.y, -4.281337229063151e-5.rad());
     }
+
+    #[test]
+    fn test_cip_celestial_to_intermediate_matrix() {
+        let xy = CipCoords {
+            x: 5.791_308_486_706_011e-4.rad(),
+            y: 4.020_579_816_732_961e-5.rad(),
+        };
+        let s = CioLocator(-1.220_040_848_472_272e-8.rad());
+        let exp = DMat3::from_cols_array(&[
+            0.999_999_832_303_715_7,
+            5.581_984_869_168_499e-10,
+            -5.791_308_491_611_282e-4,
+            -2.384_261_642_670_440_2e-8,
+            0.999_999_999_191_746_9,
+            -4.020_579_110_169_669e-5,
+            5.791_308_486_706_011e-4,
+            4.020_579_816_732_961e-5,
+            0.999_999_831_495_462_8,
+        ])
+        .transpose();
+        let act = xy.celestial_to_intermediate_matrix(s);
+        assert_approx_eq!(act, exp, atol <= 1e-12);
+    }
 }

crates/lox-earth/src/cip.rs

@@ -174,6 +174,7 @@ fn calculate_cip_unit_vector(
 #[cfg(test)]
 mod tests {
     use lox_test_utils::assert_approx_eq;
+    use lox_time::{Time, time_scales::Tdb};
 
     use crate::cip::CipCoords;
 
@@ -182,27 +183,11 @@ mod tests {
     const TOLERANCE: f64 = 1e-12;
 
     #[test]
-    fn test_cip_xy_jd0() {
-        let jd0: JulianCenturies = -67.11964407939767;
-        let CipCoords { x, y } = CipCoords::new(jd0);
-        assert_approx_eq!(x, -0.4088355637476968.rad(), rtol <= TOLERANCE);
-        assert_approx_eq!(y, -0.38359667445777073.rad(), rtol <= TOLERANCE);
-    }
-
-    #[test]
-    fn test_cip_xy_j2000() {
-        let j2000: JulianCenturies = 0.0;
-        let CipCoords { x, y } = CipCoords::new(j2000);
-        assert_approx_eq!(x, -0.0000269463795685740.rad(), rtol <= TOLERANCE);
-        assert_approx_eq!(y, -0.00002800472282281282.rad(), rtol <= TOLERANCE);
-    }
-
-    #[test]
-    fn test_cip_xy_j2100() {
-        let j2100: JulianCenturies = 1.0;
-        let CipCoords { x, y } = CipCoords::new(j2100);
-        assert_approx_eq!(x, 0.00972070446172924.rad(), rtol <= TOLERANCE);
-        assert_approx_eq!(y, -0.0000673058699616719.rad(), rtol <= TOLERANCE);
+    fn test_cip_coordinates_iau2006() {
+        let time = Time::from_two_part_julian_date(Tdb, 2400000.5, 53736.0);
+        let CipCoords { x, y } = CipCoords::iau2006(time);
+        assert_approx_eq!(x, 5.791_308_486_706_011e-4.rad(), rtol <= TOLERANCE);
+        assert_approx_eq!(y, 4.020_579_816_732_958e-5.rad(), rtol <= TOLERANCE);
     }
 
     #[test]