feat(lox-earth): add sidereal time calculations

Author: helgee

crates/lox-earth/src/itrf.rs

@@ -8,7 +8,6 @@ use crate::coordinate_transformations::{
     PoleCoords, celestial_to_intermediate_frame_of_date_matrix, polar_motion_matrix,
 };
 use crate::eop::{EopProvider, EopProviderError};
-use crate::rotation::earth_rotation_angle_00;
 use crate::tio::tio_locator;
 use glam::{DMat3, DVec3};
 use lox_bodies::{Earth, RotationalElements};
@@ -52,11 +51,12 @@ where
 }
 
 pub fn cirf_to_tirf(seconds: f64) -> Rotation {
-    let era = earth_rotation_angle_00(seconds / SECONDS_PER_DAY);
-    let rate = Earth.rotation_rate(seconds);
-    let m = DMat3::from_rotation_z(-era.to_radians());
-    let v = DVec3::new(0.0, 0.0, rate);
-    Rotation::new(m).with_angular_velocity(v)
+    // let era = earth_rotation_angle_00(seconds / SECONDS_PER_DAY);
+    // let rate = Earth.rotation_rate(seconds);
+    // let m = DMat3::from_rotation_z(-era.to_radians());
+    // let v = DVec3::new(0.0, 0.0, rate);
+    // Rotation::new(m).with_angular_velocity(v)
+    todo!()
 }
 
 pub fn tirf_to_itrf(centuries: f64) -> Rotation {

crates/lox-earth/src/lib.rs

@@ -13,6 +13,7 @@ pub mod ephemeris;
 pub mod fundamental;
 pub mod itrf;
 pub mod nutation;
+pub mod precession;
 pub mod rotation;
 #[allow(dead_code)]
 pub mod tides;

crates/lox-earth/src/precession.rs

@@ -0,0 +1,41 @@
+// SPDX-FileCopyrightText: 2025 Helge Eichhorn <[email protected]>
+//
+// SPDX-License-Identifier: MPL-2.0
+
+use lox_core::units::Angle;
+use lox_time::{Time, julian_dates::JulianDate, time_scales::Tt};
+
+use crate::nutation::Nutation;
+
+const PRECOR: Angle = Angle::arcseconds(-0.29965);
+const OBLCOR: Angle = Angle::arcseconds(-0.02524);
+
+pub type PrecessionCorrections = Nutation;
+
+pub fn precession_rate_iau2000(time: Time<Tt>) -> PrecessionCorrections {
+    let t = time.centuries_since_j2000();
+
+    PrecessionCorrections {
+        longitude: t * PRECOR,
+        obliquity: t * OBLCOR,
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use lox_core::units::AngleUnits;
+    use lox_test_utils::assert_approx_eq;
+
+    use super::*;
+
+    #[test]
+    fn test_precession_rate_iau2000() {
+        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
+        let exp = PrecessionCorrections {
+            longitude: -8.716_465_172_668_348e-8.rad(),
+            obliquity: -7.342_018_386_722_813e-9.rad(),
+        };
+        let act = precession_rate_iau2000(time);
+        assert_approx_eq!(act, exp, atol <= 1e-22);
+    }
+}

crates/lox-earth/src/rotation.rs

@@ -6,38 +6,243 @@
 //! Module rotation_angle exposes functions for calculating the Earth Rotation Angle (ERA).
 
 use std::f64::consts::TAU;
+use std::iter::zip;
 
-use lox_core::types::units::Days;
+use fast_polynomial::poly_array;
 use lox_core::units::{Angle, AngleUnits};
-use lox_time::{Time, julian_dates::JulianDate, time_scales::TimeScale};
+use lox_test_utils::ApproxEq;
+use lox_time::time_scales::{Tdb, Tt, Ut1};
+use lox_time::{Time, julian_dates::JulianDate};
 
-pub trait EarthRotation: JulianDate {
-    fn earth_rotation_angle_iau2000(&self) -> Angle;
-    fn equation_of_the_equinoxes(&self);
+use crate::ecliptic::MeanObliquity;
+use crate::fundamental::iers03::{
+    d_iers03, earth_l_iers03, f_iers03, l_iers03, lp_iers03, omega_iers03, pa_iers03,
+    venus_l_iers03,
+};
+use crate::nutation::Nutation;
+use crate::precession::precession_rate_iau2000;
+
+mod complementary_terms;
+
+#[derive(Debug, Clone, Copy, PartialOrd, PartialEq, ApproxEq)]
+pub struct EarthRotationAngle(pub Angle);
+
+impl EarthRotationAngle {
+    pub fn iau2000(time: Time<Ut1>) -> Self {
+        let d = time.days_since_j2000();
+        let f = d.rem_euclid(1.0); // fractional part of t
+        Self(
+            (TAU * (f + 0.7790572732640 + 0.00273781191135448 * d))
+                .rad()
+                .mod_two_pi(),
+        )
+    }
+}
+
+pub type Era = EarthRotationAngle;
+
+#[derive(Debug, Clone, Copy, PartialOrd, PartialEq, ApproxEq)]
+pub struct GreenwichMeanSiderealTime(pub Angle);
+
+pub type Gmst = GreenwichMeanSiderealTime;
+
+impl GreenwichMeanSiderealTime {
+    pub fn iau2000(tt: Time<Tt>, ut1: Time<Ut1>) -> Self {
+        let t = tt.centuries_since_j2000();
+        Self(
+            EarthRotationAngle::iau2000(ut1).0
+                + Angle::arcseconds(poly_array(
+                    t,
+                    &[0.014506, 4612.15739966, 1.39667721, -0.00009344, 0.00001882],
+                ))
+                .mod_two_pi(),
+        )
+    }
+
+    pub fn iau2006(tt: Time<Tt>, ut1: Time<Ut1>) -> Self {
+        let t = tt.centuries_since_j2000();
+        Self(
+            EarthRotationAngle::iau2000(ut1).0
+                + Angle::arcseconds(poly_array(
+                    t,
+                    &[
+                        0.014506,
+                        4612.156534,
+                        1.3915817,
+                        -0.00000044,
+                        -0.000029956,
+                        -0.0000000368,
+                    ],
+                ))
+                .mod_two_pi(),
+        )
+    }
 }
 
-// impl<T> EarthRotation for Time<T> where T: TimeScale {}
+#[derive(Debug, Clone, Copy, PartialOrd, PartialEq, ApproxEq)]
+pub struct GreenwichApparentSiderealTime(pub Angle);
 
-pub fn earth_rotation_angle_00(days_since_j2000_ut1: Days) -> Angle {
-    let f = days_since_j2000_ut1.rem_euclid(1.0); // fractional part of t
-    let era = (TAU * (f + 0.7790572732640 + 0.00273781191135448 * days_since_j2000_ut1)).rad();
-    era.mod_two_pi()
+pub type Gast = GreenwichApparentSiderealTime;
+
+impl GreenwichApparentSiderealTime {
+    pub fn iau2000a(tt: Time<Tt>, ut1: Time<Ut1>) -> Self {
+        Self(
+            (GreenwichMeanSiderealTime::iau2000(tt, ut1).0
+                + EquationOfTheEquinoxes::iau2000a(tt).0)
+                .mod_two_pi(),
+        )
+    }
+
+    pub fn iau2000b(tt: Time<Tt>, ut1: Time<Ut1>) -> Self {
+        Self(
+            (GreenwichMeanSiderealTime::iau2000(tt, ut1).0
+                + EquationOfTheEquinoxes::iau2000b(tt).0)
+                .mod_two_pi(),
+        )
+    }
+}
+
+#[derive(Debug, Clone, Copy, PartialOrd, PartialEq, ApproxEq)]
+pub struct EquationOfTheEquinoxes(pub Angle);
+
+impl EquationOfTheEquinoxes {
+    pub fn iau2000a(time: Time<Tt>) -> Self {
+        let Nutation {
+            obliquity: depspr, ..
+        } = precession_rate_iau2000(time);
+        let epsa = MeanObliquity::iau1980(time).0 + depspr;
+        let Nutation {
+            longitude: dpsi, ..
+        } = Nutation::iau2000a(time.with_scale(Tdb));
+        Self::iau2000(time, epsa, dpsi)
+    }
+
+    pub fn iau2000b(time: Time<Tt>) -> Self {
+        let Nutation {
+            obliquity: depspr, ..
+        } = precession_rate_iau2000(time);
+        let epsa = MeanObliquity::iau1980(time).0 + depspr;
+        let Nutation {
+            longitude: dpsi, ..
+        } = Nutation::iau2000b(time.with_scale(Tdb));
+        Self::iau2000(time, epsa, dpsi)
+    }
+
+    pub fn iau2000(time: Time<Tt>, epsa: Angle, dpsi: Angle) -> Self {
+        Self(epsa.cos() * dpsi + Self::complimentary_terms_iau2000(time))
+    }
+
+    fn complimentary_terms_iau2000(time: Time<Tt>) -> Angle {
+        let t = time.centuries_since_j2000();
+
+        let fa = [
+            l_iers03(t).as_f64(),
+            lp_iers03(t).as_f64(),
+            f_iers03(t).as_f64(),
+            d_iers03(t).as_f64(),
+            omega_iers03(t).as_f64(),
+            venus_l_iers03(t).as_f64(),
+            earth_l_iers03(t).as_f64(),
+            pa_iers03(t).as_f64(),
+        ];
+
+        let s0 = complementary_terms::E0.iter().rev().fold(0.0, |s0, term| {
+            let a = zip(&term.nfa, &fa).fold(0.0, |a, (&nfa, &fa)| a + nfa as f64 * fa);
+            let (sa, ca) = a.sin_cos();
+            s0 + term.s * sa + term.c * ca
+        });
+
+        let s1 = complementary_terms::E1.iter().rev().fold(0.0, |s1, term| {
+            let a = zip(&term.nfa, &fa).fold(0.0, |a, (&nfa, &fa)| a + nfa as f64 * fa);
+            let (sa, ca) = a.sin_cos();
+            s1 + term.s * sa + term.c * ca
+        });
+
+        Angle::arcseconds(s0 + s1 * t)
+    }
 }
 
 #[cfg(test)]
 mod tests {
     use lox_test_utils::assert_approx_eq;
 
-    use rstest::rstest;
-
     use super::*;
 
-    #[rstest]
-    #[case::before_j2000(-123.45, 6.227104062035152.rad())]
-    #[case::j2000(0.0, 4.894961212823756.rad())]
-    #[case::after_j2000(123.45, 3.562818363612361.rad())]
-    fn test_rotation_angle_00(#[case] days_since_j2000_ut1: Days, #[case] expected: Angle) {
-        let actual = earth_rotation_angle_00(days_since_j2000_ut1);
-        assert_approx_eq!(expected, actual, rtol <= 1e-9);
+    #[test]
+    fn test_earth_rotation_angle_iau2000() {
+        let time = Time::from_two_part_julian_date(Ut1, 2400000.5, 54388.0);
+        let exp = EarthRotationAngle(Angle::new(0.402_283_724_002_815_8));
+        let act = EarthRotationAngle::iau2000(time);
+        assert_approx_eq!(act, exp, rtol <= 1e-12);
+    }
+
+    #[test]
+    fn test_greenwich_apparent_sidereal_time_iau2000a() {
+        let tt = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
+        let ut1 = Time::from_two_part_julian_date(Ut1, 2400000.5, 53736.0);
+        let exp = GreenwichApparentSiderealTime(Angle::new(1.754_166_138_018_281_4));
+        let act = Gast::iau2000a(tt, ut1);
+        assert_approx_eq!(act, exp, rtol <= 1e-12);
+    }
+
+    #[test]
+    fn test_greenwich_apparent_sidereal_time_iau2000b() {
+        let tt = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
+        let ut1 = Time::from_two_part_julian_date(Ut1, 2400000.5, 53736.0);
+        let exp = GreenwichApparentSiderealTime(Angle::new(1.754_166_136_510_680_7));
+        let act = Gast::iau2000b(tt, ut1);
+        assert_approx_eq!(act, exp, rtol <= 1e-12);
+    }
+
+    #[test]
+    fn test_greenwich_mean_sidereal_time_iau2000() {
+        let tt = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
+        let ut1 = Time::from_two_part_julian_date(Ut1, 2400000.5, 53736.0);
+        let exp = GreenwichMeanSiderealTime(Angle::new(1.754_174_972_210_740_7));
+        let act = Gmst::iau2000(tt, ut1);
+        assert_approx_eq!(act, exp, rtol <= 1e-12);
+    }
+
+    #[test]
+    fn test_greenwich_mean_sidereal_time_iau2006() {
+        let tt = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
+        let ut1 = Time::from_two_part_julian_date(Ut1, 2400000.5, 53736.0);
+        let exp = GreenwichMeanSiderealTime(Angle::new(1.754_174_971_870_091_2));
+        let act = Gmst::iau2006(tt, ut1);
+        assert_approx_eq!(act, exp, rtol <= 1e-12);
+    }
+
+    #[test]
+    fn test_equation_of_the_equinoxes_iau2000a() {
+        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
+        let exp = EquationOfTheEquinoxes(Angle::new(-8.834_192_459_222_587e-6));
+        let act = EquationOfTheEquinoxes::iau2000a(time);
+        assert_approx_eq!(act, exp, atol <= 1e-18);
+    }
+
+    #[test]
+    fn test_equation_of_the_equinoxes_iau2000b() {
+        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
+        let exp = EquationOfTheEquinoxes(Angle::new(-8.835_700_060_003_032e-6));
+        let act = EquationOfTheEquinoxes::iau2000b(time);
+        assert_approx_eq!(act, exp, atol <= 1e-18);
+    }
+
+    #[test]
+    fn test_equation_of_the_equinoxes_iau2000() {
+        let epsa = Angle::new(0.409_078_976_335_651);
+        let dpsi = Angle::new(-9.630_909_107_115_582e-6);
+        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
+        let exp = EquationOfTheEquinoxes(Angle::new(-8.834_193_235_367_966e-6));
+        let act = EquationOfTheEquinoxes::iau2000(time, epsa, dpsi);
+        assert_approx_eq!(act, exp, atol <= 1e-20);
+    }
+
+    #[test]
+    fn test_equation_of_the_equinoxes_complimentary_terms() {
+        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
+        let exp = Angle::new(2.046_085_004_885_125e-9);
+        let act = EquationOfTheEquinoxes::complimentary_terms_iau2000(time);
+        assert_approx_eq!(act, exp, atol <= 1e-20);
     }
 }