lox_frames/iers/
earth_rotation.rs

1// SPDX-FileCopyrightText: 2024 Angus Morrison <github@angus-morrison.com>
2// SPDX-FileCopyrightText: 2024 Helge Eichhorn <git@helgeeichhorn.de>
3//
4// SPDX-License-Identifier: MPL-2.0
5
6//! Module rotation_angle exposes functions for calculating the Earth Rotation Angle (ERA).
7
8use std::f64::consts::TAU;
9use std::iter::zip;
10
11use fast_polynomial::poly_array;
12use lox_approx::ApproxEq;
13use lox_core::f64::consts::{SECONDS_PER_DAY, SECONDS_PER_HALF_DAY};
14use lox_core::glam::DMat3;
15use lox_core::units::{Angle, AngleUnits, Sign};
16use lox_time::time_scales::{Tdb, Tt, Ut1};
17use lox_time::{Time, julian_dates::JulianDate};
18
19use crate::iers::ecliptic::MeanObliquity;
20use crate::iers::fundamental::iers03::{
21    d_iers03, earth_l_iers03, f_iers03, l_iers03, lp_iers03, omega_iers03, pa_iers03,
22    venus_l_iers03,
23};
24use crate::iers::nutation::Nutation;
25use crate::iers::precession::PrecessionCorrectionsIau2000;
26use crate::iers::{Corrections, Iau2000Model, ReferenceSystem};
27
28mod complementary_terms;
29
30impl ReferenceSystem {
31    /// Returns the Earth rotation matrix for the given times and corrections.
32    pub fn earth_rotation(&self, tt: Time<Tt>, ut1: Time<Ut1>, corr: Corrections) -> DMat3 {
33        self.greenwich_apparent_sidereal_time(tt, ut1, corr)
34            .0
35            .rotation_z()
36    }
37
38    /// Computes Greenwich Mean Sidereal Time (GMST) for the active IERS convention.
39    pub fn greenwich_mean_sidereal_time(
40        &self,
41        tt: Time<Tt>,
42        ut1: Time<Ut1>,
43    ) -> GreenwichMeanSiderealTime {
44        match self {
45            ReferenceSystem::Iers1996 => GreenwichMeanSiderealTime::iau1982(ut1),
46            ReferenceSystem::Iers2003(_) => GreenwichMeanSiderealTime::iau2000(tt, ut1),
47            ReferenceSystem::Iers2010 => GreenwichMeanSiderealTime::iau2006(tt, ut1),
48        }
49    }
50
51    /// Computes Greenwich Apparent Sidereal Time (GAST) for the active IERS convention.
52    pub fn greenwich_apparent_sidereal_time(
53        &self,
54        tt: Time<Tt>,
55        ut1: Time<Ut1>,
56        corr: Corrections,
57    ) -> GreenwichApparentSiderealTime {
58        if corr.is_zero() {
59            return match self {
60                ReferenceSystem::Iers1996 => Gast::iau1994(ut1),
61                ReferenceSystem::Iers2003(model) => match model {
62                    Iau2000Model::A => Gast::iau2000a(tt, ut1),
63                    Iau2000Model::B => Gast::iau2000b(tt, ut1),
64                },
65                ReferenceSystem::Iers2010 => Gast::iau2006a(tt, ut1),
66            };
67        };
68        let gmst = self.greenwich_mean_sidereal_time(tt, ut1);
69        let tdb = tt.with_scale(Tdb);
70        let rpb = self.bias_precession_matrix(tt);
71        let epsa = self.mean_obliquity(tt);
72        let mut nut = self.nutation(tdb);
73        let ecl_corr = self.ecliptic_corrections(corr, nut, epsa, rpb);
74        nut += ecl_corr;
75        match self {
76            ReferenceSystem::Iers1996 => {
77                let ee = EquationOfTheEquinoxes::iau1994(tdb);
78                GreenwichApparentSiderealTime(
79                    (gmst.0 + ee.0 + epsa.0.cos() * ecl_corr.0).mod_two_pi(),
80                )
81            }
82            ReferenceSystem::Iers2003(_) => {
83                let ee = EquationOfTheEquinoxes::iau2000(tt, epsa.0, nut.dpsi);
84                GreenwichApparentSiderealTime((gmst.0 + ee.0).mod_two_pi())
85            }
86            ReferenceSystem::Iers2010 => {
87                let ee = EquationOfTheEquinoxes::iau2000(tt, epsa.0, nut.dpsi);
88                GreenwichApparentSiderealTime((gmst.0 + ee.0).mod_two_pi())
89            }
90        }
91    }
92}
93
94/// Earth Rotation Angle (ERA).
95#[derive(Debug, Clone, Copy, PartialOrd, PartialEq, ApproxEq)]
96pub struct EarthRotationAngle(pub Angle);
97
98impl EarthRotationAngle {
99    /// Computes the ERA using the IAU 2000 model.
100    pub fn iau2000(time: Time<Ut1>) -> Self {
101        let d = time.days_since_j2000();
102        let f = d.rem_euclid(1.0); // fractional part of t
103        Self(
104            (TAU * (f + 0.7790572732640 + 0.00273781191135448 * d))
105                .rad()
106                .mod_two_pi(),
107        )
108    }
109}
110
111/// Alias for [`EarthRotationAngle`].
112pub type Era = EarthRotationAngle;
113
114/// Greenwich Mean Sidereal Time (GMST).
115#[derive(Debug, Clone, Copy, PartialOrd, PartialEq, ApproxEq)]
116pub struct GreenwichMeanSiderealTime(pub Angle);
117
118/// Alias for [`GreenwichMeanSiderealTime`].
119pub type Gmst = GreenwichMeanSiderealTime;
120
121// Coefficients of IAU 1982 GMST-UT1 model
122const A: f64 = 24110.54841 - SECONDS_PER_HALF_DAY;
123const B: f64 = 8640184.812866;
124const C: f64 = 0.093104;
125const D: f64 = -6.2e-6;
126
127impl GreenwichMeanSiderealTime {
128    /// Computes GMST using the IAU 1982 model.
129    pub fn iau1982(time: Time<Ut1>) -> Self {
130        let t = time.centuries_since_j2000();
131        let f = time.days_since_j2000().rem_euclid(1.0) * SECONDS_PER_DAY;
132        Self(Angle::from_hms(Sign::Positive, 0, 0, poly_array(t, &[A, B, C, D]) + f).mod_two_pi())
133    }
134
135    /// Computes GMST using the IAU 2000 model.
136    pub fn iau2000(tt: Time<Tt>, ut1: Time<Ut1>) -> Self {
137        let t = tt.centuries_since_j2000();
138        Self(
139            EarthRotationAngle::iau2000(ut1).0
140                + Angle::arcseconds(poly_array(
141                    t,
142                    &[0.014506, 4612.15739966, 1.39667721, -0.00009344, 0.00001882],
143                ))
144                .mod_two_pi(),
145        )
146    }
147
148    /// Computes GMST using the IAU 2006 model.
149    pub fn iau2006(tt: Time<Tt>, ut1: Time<Ut1>) -> Self {
150        let t = tt.centuries_since_j2000();
151        Self(
152            EarthRotationAngle::iau2000(ut1).0
153                + Angle::arcseconds(poly_array(
154                    t,
155                    &[
156                        0.014506,
157                        4612.156534,
158                        1.3915817,
159                        -0.00000044,
160                        -0.000029956,
161                        -0.0000000368,
162                    ],
163                ))
164                .mod_two_pi(),
165        )
166    }
167}
168
169/// Greenwich Apparent Sidereal Time (GAST).
170#[derive(Debug, Clone, Copy, PartialOrd, PartialEq, ApproxEq)]
171pub struct GreenwichApparentSiderealTime(pub Angle);
172
173/// Alias for [`GreenwichApparentSiderealTime`].
174pub type Gast = GreenwichApparentSiderealTime;
175
176impl GreenwichApparentSiderealTime {
177    /// Computes GAST using the IAU 1994 model.
178    pub fn iau1994(time: Time<Ut1>) -> Self {
179        Self(
180            (GreenwichMeanSiderealTime::iau1982(time).0
181                + EquationOfTheEquinoxes::iau1994(time.with_scale(Tdb)).0)
182                .mod_two_pi(),
183        )
184    }
185
186    /// Computes GAST using the IAU 2000A model.
187    pub fn iau2000a(tt: Time<Tt>, ut1: Time<Ut1>) -> Self {
188        Self(
189            (GreenwichMeanSiderealTime::iau2000(tt, ut1).0
190                + EquationOfTheEquinoxes::iau2000a(tt).0)
191                .mod_two_pi(),
192        )
193    }
194
195    /// Computes GAST using the IAU 2000B model.
196    pub fn iau2000b(tt: Time<Tt>, ut1: Time<Ut1>) -> Self {
197        Self(
198            (GreenwichMeanSiderealTime::iau2000(tt, ut1).0
199                + EquationOfTheEquinoxes::iau2000b(tt).0)
200                .mod_two_pi(),
201        )
202    }
203
204    /// Computes GAST using the IAU 2006A model.
205    pub fn iau2006a(tt: Time<Tt>, ut1: Time<Ut1>) -> Self {
206        Self(
207            (GreenwichMeanSiderealTime::iau2006(tt, ut1).0
208                + EquationOfTheEquinoxes::iau2006a(tt).0)
209                .mod_two_pi(),
210        )
211    }
212}
213
214/// Equation of the equinoxes (GAST - GMST).
215#[derive(Debug, Clone, Copy, PartialOrd, PartialEq, ApproxEq)]
216pub struct EquationOfTheEquinoxes(pub Angle);
217
218impl EquationOfTheEquinoxes {
219    /// Computes the equation of the equinoxes using the IAU 1994 model.
220    pub fn iau1994(time: Time<Tdb>) -> Self {
221        let Nutation { dpsi, .. } = Nutation::iau1980(time);
222        Self::iau1994_from_dpsi(time, dpsi)
223    }
224
225    /// Computes the IAU 1994 equation of the equinoxes from a precomputed
226    /// nutation in longitude.
227    pub(crate) fn iau1994_from_dpsi(time: Time<Tdb>, dpsi: Angle) -> Self {
228        let t = time.centuries_since_j2000();
229        let om = (Angle::arcseconds(poly_array(t, &[450160.280, -482890.539, 7.455, 0.008]))
230            + Angle::new((-5.0 * t).rem_euclid(1.0) * TAU))
231        .mod_two_pi();
232        let eps0 = MeanObliquity::iau1980(time.with_scale(Tt));
233        Self(
234            eps0.0.cos() * dpsi
235                + Angle::arcseconds(0.00264 * om.sin() + 0.000063 * (om + om).sin()),
236        )
237    }
238
239    /// Computes the equation of the equinoxes using the IAU 2000A model.
240    pub fn iau2000a(time: Time<Tt>) -> Self {
241        let PrecessionCorrectionsIau2000 { depspr, .. } = PrecessionCorrectionsIau2000::new(time);
242        let epsa = MeanObliquity::iau1980(time).0 + depspr;
243        let Nutation { dpsi, .. } = Nutation::iau2000a(time.with_scale(Tdb));
244        Self::iau2000(time, epsa, dpsi)
245    }
246
247    /// Computes the equation of the equinoxes using the IAU 2000B model.
248    pub fn iau2000b(time: Time<Tt>) -> Self {
249        let PrecessionCorrectionsIau2000 { depspr, .. } = PrecessionCorrectionsIau2000::new(time);
250        let epsa = MeanObliquity::iau1980(time).0 + depspr;
251        let Nutation { dpsi, .. } = Nutation::iau2000b(time.with_scale(Tdb));
252        Self::iau2000(time, epsa, dpsi)
253    }
254
255    /// Computes the equation of the equinoxes using the IAU 2006A model.
256    pub fn iau2006a(time: Time<Tt>) -> Self {
257        let epsa = MeanObliquity::iau2006(time);
258        let Nutation { dpsi, .. } = Nutation::iau2006a(time.with_scale(Tdb));
259        Self::iau2000(time, epsa.0, dpsi)
260    }
261
262    /// Computes the equation of the equinoxes given mean obliquity and nutation in longitude.
263    pub fn iau2000(time: Time<Tt>, epsa: Angle, dpsi: Angle) -> Self {
264        Self(epsa.cos() * dpsi + Self::complimentary_terms_iau2000(time))
265    }
266
267    fn complimentary_terms_iau2000(time: Time<Tt>) -> Angle {
268        let t = time.centuries_since_j2000();
269
270        let fa = [
271            l_iers03(t).as_f64(),
272            lp_iers03(t).as_f64(),
273            f_iers03(t).as_f64(),
274            d_iers03(t).as_f64(),
275            omega_iers03(t).as_f64(),
276            venus_l_iers03(t).as_f64(),
277            earth_l_iers03(t).as_f64(),
278            pa_iers03(t).as_f64(),
279        ];
280
281        let s0 = complementary_terms::E0.iter().rev().fold(0.0, |s0, term| {
282            let a = zip(&term.nfa, &fa).fold(0.0, |a, (&nfa, &fa)| a + nfa as f64 * fa);
283            let (sa, ca) = a.sin_cos();
284            s0 + term.s * sa + term.c * ca
285        });
286
287        let s1 = complementary_terms::E1.iter().rev().fold(0.0, |s1, term| {
288            let a = zip(&term.nfa, &fa).fold(0.0, |a, (&nfa, &fa)| a + nfa as f64 * fa);
289            let (sa, ca) = a.sin_cos();
290            s1 + term.s * sa + term.c * ca
291        });
292
293        Angle::arcseconds(s0 + s1 * t)
294    }
295}
296
297#[cfg(test)]
298mod tests {
299    use lox_approx::assert_approx_eq;
300
301    use super::*;
302
303    #[test]
304    fn test_earth_rotation_angle_iau2000() {
305        let time = Time::from_two_part_julian_date(Ut1, 2400000.5, 54388.0);
306        let exp = EarthRotationAngle(Angle::new(0.402_283_724_002_815_8));
307        let act = EarthRotationAngle::iau2000(time);
308        assert_approx_eq!(act, exp, rtol <= 1e-12);
309    }
310
311    #[test]
312    fn test_greenwich_apparent_sidereal_time_iau1994() {
313        let ut1 = Time::from_two_part_julian_date(Ut1, 2400000.5, 53736.0);
314        let exp = GreenwichApparentSiderealTime(Angle::new(1.754_166_136_020_645_3));
315        let act = Gast::iau1994(ut1);
316        assert_approx_eq!(act, exp, rtol <= 1e-12);
317    }
318
319    #[test]
320    fn test_greenwich_apparent_sidereal_time_iau2000a() {
321        let tt = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
322        let ut1 = Time::from_two_part_julian_date(Ut1, 2400000.5, 53736.0);
323        let exp = GreenwichApparentSiderealTime(Angle::new(1.754_166_138_018_281_4));
324        let act = Gast::iau2000a(tt, ut1);
325        assert_approx_eq!(act, exp, rtol <= 1e-12);
326    }
327
328    #[test]
329    fn test_greenwich_apparent_sidereal_time_iau2000b() {
330        let tt = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
331        let ut1 = Time::from_two_part_julian_date(Ut1, 2400000.5, 53736.0);
332        let exp = GreenwichApparentSiderealTime(Angle::new(1.754_166_136_510_680_7));
333        let act = Gast::iau2000b(tt, ut1);
334        assert_approx_eq!(act, exp, rtol <= 1e-12);
335    }
336
337    #[test]
338    fn test_greenwich_mean_sidereal_time_iau1982() {
339        let ut1 = Time::from_two_part_julian_date(Ut1, 2400000.5, 53736.0);
340        let exp = GreenwichMeanSiderealTime(Angle::new(1.754_174_981_860_675));
341        let act = Gmst::iau1982(ut1);
342        assert_approx_eq!(act, exp, rtol <= 1e-12);
343    }
344
345    #[test]
346    fn test_greenwich_mean_sidereal_time_iau2000() {
347        let tt = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
348        let ut1 = Time::from_two_part_julian_date(Ut1, 2400000.5, 53736.0);
349        let exp = GreenwichMeanSiderealTime(Angle::new(1.754_174_972_210_740_7));
350        let act = Gmst::iau2000(tt, ut1);
351        assert_approx_eq!(act, exp, rtol <= 1e-12);
352    }
353
354    #[test]
355    fn test_greenwich_mean_sidereal_time_iau2006() {
356        let tt = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
357        let ut1 = Time::from_two_part_julian_date(Ut1, 2400000.5, 53736.0);
358        let exp = GreenwichMeanSiderealTime(Angle::new(1.754_174_971_870_091_2));
359        let act = Gmst::iau2006(tt, ut1);
360        assert_approx_eq!(act, exp, rtol <= 1e-12);
361    }
362
363    #[test]
364    fn test_equation_of_the_equinoxes_iau1994() {
365        let time = Time::from_two_part_julian_date(Tdb, 2400000.5, 41234.0);
366        let exp = EquationOfTheEquinoxes(Angle::new(5.357_758_254_609_257e-5));
367        let act = EquationOfTheEquinoxes::iau1994(time);
368        assert_approx_eq!(act, exp, atol <= 1e-17);
369    }
370
371    #[test]
372    fn test_equation_of_the_equinoxes_iau2000a() {
373        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
374        let exp = EquationOfTheEquinoxes(Angle::new(-8.834_192_459_222_587e-6));
375        let act = EquationOfTheEquinoxes::iau2000a(time);
376        assert_approx_eq!(act, exp, atol <= 1e-18);
377    }
378
379    #[test]
380    fn test_equation_of_the_equinoxes_iau2000b() {
381        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
382        let exp = EquationOfTheEquinoxes(Angle::new(-8.835_700_060_003_032e-6));
383        let act = EquationOfTheEquinoxes::iau2000b(time);
384        assert_approx_eq!(act, exp, atol <= 1e-18);
385    }
386
387    #[test]
388    fn test_equation_of_the_equinoxes_iau2000() {
389        let epsa = Angle::new(0.409_078_976_335_651);
390        let dpsi = Angle::new(-9.630_909_107_115_582e-6);
391        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
392        let exp = EquationOfTheEquinoxes(Angle::new(-8.834_193_235_367_966e-6));
393        let act = EquationOfTheEquinoxes::iau2000(time, epsa, dpsi);
394        assert_approx_eq!(act, exp, atol <= 1e-20);
395    }
396
397    #[test]
398    fn test_equation_of_the_equinoxes_complimentary_terms() {
399        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
400        let exp = Angle::new(2.046_085_004_885_125e-9);
401        let act = EquationOfTheEquinoxes::complimentary_terms_iau2000(time);
402        assert_approx_eq!(act, exp, atol <= 1e-20);
403    }
404
405    #[test]
406    fn test_greenwich_apparent_sidereal_time_iau2006a() {
407        let tt = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
408        let ut1 = Time::from_two_part_julian_date(Ut1, 2400000.5, 53736.0);
409        let exp = GreenwichApparentSiderealTime(Angle::new(1.754_166_137_675_019_2));
410        let act = Gast::iau2006a(tt, ut1);
411        assert_approx_eq!(act, exp, rtol <= 1e-12);
412    }
413
414    #[test]
415    fn test_equation_of_the_equinoxes_iau2006a() {
416        // ERFA ee06a reference value. ERFA computes this via gst06a - gmst06 which
417        // takes a different computational path (through the full NPB matrix) than our
418        // direct epsa.cos() * dpsi + complementary_terms. The ~8e-13 rad difference
419        // (~0.2 μas) is expected.
420        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
421        let exp = EquationOfTheEquinoxes(Angle::new(-8.834_195_072_043_79e-6));
422        let act = EquationOfTheEquinoxes::iau2006a(time);
423        assert_approx_eq!(act, exp, atol <= 1e-12);
424    }
425}
lox_frames/iers/earth_rotation.rs

lox_frames/iers/
earth_rotation.rs
