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

lox_frames/iers/
earth_rotation.rs
