lox_frames/iers/
precession.rs

1// SPDX-FileCopyrightText: 2025 Helge Eichhorn <git@helgeeichhorn.de>
2//
3// SPDX-License-Identifier: MPL-2.0
4
5use fast_polynomial::poly_array;
6use glam::DMat3;
7use lox_test_utils::ApproxEq;
8use lox_time::{Time, julian_dates::JulianDate, time_scales::Tt};
9use lox_units::Angle;
10
11use crate::iers::{ReferenceSystem, ecliptic::MeanObliquity};
12
13impl ReferenceSystem {
14    pub fn precession(&self, time: Time<Tt>) -> PrecessionAngles {
15        match self {
16            ReferenceSystem::Iers1996 => PrecessionAngles::iau1976(time),
17            ReferenceSystem::Iers2003(_) => PrecessionAngles::iau2000(time),
18            ReferenceSystem::Iers2010 => PrecessionAngles::iau2006(time),
19        }
20    }
21    /// Returns the rotation from ICRF to the Mean-of-Date (MOD) frame for the given IERS convention.
22    pub fn bias_precession_matrix(&self, time: Time<Tt>) -> DMat3 {
23        self.precession(time).bias_precession_matrix()
24    }
25
26    /// Returns the precession-only rotation from J2000 to the Mean-of-Date (MOD) frame.
27    pub fn precession_matrix(&self, time: Time<Tt>) -> DMat3 {
28        self.precession(time).precession_matrix()
29    }
30}
31
32pub enum PrecessionAngles {
33    /// IAU1976 precession angles.
34    Iau1976 { zeta: Angle, z: Angle, theta: Angle },
35    /// IAU2000 precession angles.
36    Iau2000 {
37        psia: Angle,
38        oma: Angle,
39        chia: Angle,
40    },
41    /// IAU2006 Fukushima-Williams precession angles.
42    Iau2006 {
43        gamb: Angle,
44        phib: Angle,
45        psib: Angle,
46        epsa: Angle,
47    },
48}
49
50impl PrecessionAngles {
51    pub fn iau1976(time: Time<Tt>) -> Self {
52        let t = time.centuries_since_j2000();
53
54        let tas2r = Angle::arcseconds(t);
55        let w = 2306.2181;
56        let zeta = (w + ((0.30188) + 0.017998 * t) * t) * tas2r;
57        let z = (w + ((1.09468) + 0.018203 * t) * t) * tas2r;
58        let theta = (2004.3109 + ((-0.42665) - 0.041833 * t) * t) * tas2r;
59
60        Self::Iau1976 { zeta, z, theta }
61    }
62
63    pub fn iau2000(time: Time<Tt>) -> Self {
64        let t = time.centuries_since_j2000();
65
66        // Precession angles (Lieske et al. 1977)
67        let psia77 = Angle::arcseconds(poly_array(t, &[0.0, 5038.7784, -1.07259, -0.001147]));
68        let oma77 = EPS0 + Angle::arcseconds(poly_array(t, &[0.0, 0.0, 0.05127, -0.007726]));
69        let chia = Angle::arcseconds(poly_array(t, &[0.0, 10.5526, -2.38064, -0.001125]));
70
71        // Apply IAU 2000 precession corrections
72        let PrecessionCorrectionsIau2000 { dpsipr, depspr } =
73            PrecessionCorrectionsIau2000::new(time);
74        let psia = psia77 + dpsipr;
75        let oma = oma77 + depspr;
76
77        Self::Iau2000 { psia, oma, chia }
78    }
79
80    pub fn iau2006(time: Time<Tt>) -> Self {
81        let t = time.centuries_since_j2000();
82
83        let gamb = Angle::arcseconds(poly_array(
84            t,
85            &[
86                -0.052928,
87                10.556378,
88                0.4932044,
89                -0.00031238,
90                -0.000002788,
91                0.0000000260,
92            ],
93        ));
94        let phib = Angle::arcseconds(poly_array(
95            t,
96            &[
97                84381.412819,
98                -46.811016,
99                0.0511268,
100                0.00053289,
101                -0.000000440,
102                -0.0000000176,
103            ],
104        ));
105        let psib = Angle::arcseconds(poly_array(
106            t,
107            &[
108                -0.041775,
109                5038.481484,
110                1.5584175,
111                -0.00018522,
112                -0.000026452,
113                -0.0000000148,
114            ],
115        ));
116        let epsa = MeanObliquity::iau2006(time.with_scale(Tt)).0;
117
118        Self::Iau2006 {
119            gamb,
120            phib,
121            psib,
122            epsa,
123        }
124    }
125
126    pub fn bias_precession_matrix(&self) -> DMat3 {
127        match *self {
128            PrecessionAngles::Iau1976 { zeta, z, theta } => {
129                let rp = (-z).rotation_z() * theta.rotation_y() * (-zeta).rotation_z();
130                rp * frame_bias()
131            }
132            PrecessionAngles::Iau2000 { psia, oma, chia } => {
133                let rp = chia.rotation_z()
134                    * (-oma).rotation_x()
135                    * (-psia).rotation_z()
136                    * EPS0.rotation_x();
137
138                rp * frame_bias()
139            }
140            PrecessionAngles::Iau2006 {
141                gamb,
142                phib,
143                psib,
144                epsa,
145            } => {
146                (-epsa).rotation_x() * (-psib).rotation_z() * phib.rotation_x() * gamb.rotation_z()
147            }
148        }
149    }
150
151    /// Returns the precession-only rotation from J2000 to the Mean-of-Date (MOD) frame,
152    /// without the frame bias.
153    pub fn precession_matrix(&self) -> DMat3 {
154        match *self {
155            PrecessionAngles::Iau1976 { zeta, z, theta } => {
156                (-z).rotation_z() * theta.rotation_y() * (-zeta).rotation_z()
157            }
158            PrecessionAngles::Iau2000 { psia, oma, chia } => {
159                chia.rotation_z() * (-oma).rotation_x() * (-psia).rotation_z() * EPS0.rotation_x()
160            }
161            PrecessionAngles::Iau2006 {
162                gamb,
163                phib,
164                psib,
165                epsa,
166            } => {
167                let bp = (-epsa).rotation_x()
168                    * (-psib).rotation_z()
169                    * phib.rotation_x()
170                    * gamb.rotation_z();
171                bp * frame_bias().transpose()
172            }
173        }
174    }
175}
176
177const D_PSI_BIAS: Angle = Angle::arcseconds(-0.041775);
178const D_EPS_BIAS: Angle = Angle::arcseconds(-0.0068192);
179const D_RA0: Angle = Angle::arcseconds(-0.0146);
180const EPS0: Angle = Angle::arcseconds(84381.448);
181
182pub fn frame_bias() -> DMat3 {
183    (-D_EPS_BIAS).rotation_x() * (EPS0.sin() * D_PSI_BIAS).rotation_y() * D_RA0.rotation_z()
184}
185
186const PRECOR: Angle = Angle::arcseconds(-0.29965);
187const OBLCOR: Angle = Angle::arcseconds(-0.02524);
188
189#[derive(Debug, Clone, Copy, PartialEq, PartialOrd, ApproxEq)]
190pub struct PrecessionCorrectionsIau2000 {
191    pub dpsipr: Angle,
192    pub depspr: Angle,
193}
194
195impl PrecessionCorrectionsIau2000 {
196    pub fn new(time: Time<Tt>) -> Self {
197        let t = time.centuries_since_j2000();
198
199        Self {
200            dpsipr: t * PRECOR,
201            depspr: t * OBLCOR,
202        }
203    }
204}
205
206#[cfg(test)]
207mod tests {
208    use lox_test_utils::assert_approx_eq;
209    use lox_units::AngleUnits;
210
211    use crate::iers::Iau2000Model;
212
213    use super::*;
214
215    #[test]
216    fn test_precession_iers1996() {
217        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 50123.9999);
218        let exp = DMat3::from_cols_array(&[
219            0.999_999_550_432_835_1,
220            8.696_632_209_480_961e-4,
221            3.779_153_474_959_888_4e-4,
222            -8.696_632_209_485_112e-4,
223            0.999_999_621_842_856_1,
224            -1.643_284_776_111_886_4e-7,
225            -3.779_153_474_950_335e-4,
226            -1.643_306_746_147_367e-7,
227            0.999_999_928_589_979,
228        ])
229        .transpose()
230            * frame_bias();
231        let act = ReferenceSystem::Iers1996.bias_precession_matrix(time);
232        assert_approx_eq!(act, exp, atol <= 1e-12);
233    }
234
235    #[test]
236    fn test_bias_precession_iau2000() {
237        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 50123.9999);
238        let exp = DMat3::from_cols_array(&[
239            0.999_999_550_517_508_8,
240            8.695_405_883_617_885e-4,
241            3.779_734_722_239_007e-4,
242            -8.695_405_990_410_864e-4,
243            0.999_999_621_949_492_5,
244            -1.360_775_820_404_982e-7,
245            -3.779_734_476_558_185e-4,
246            -1.925_857_585_832_024e-7,
247            0.999_999_928_568_015_3,
248        ])
249        .transpose();
250        let act = ReferenceSystem::Iers2003(Iau2000Model::A).bias_precession_matrix(time);
251        assert_approx_eq!(act, exp, atol <= 1e-12);
252    }
253
254    #[test]
255    fn test_precession_iau2006() {
256        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 50123.9999);
257        let exp = DMat3::from_cols_array(&[
258            0.999_999_550_517_600_7,
259            8.695_404_617_348_209e-4,
260            3.779_735_201_865_589e-4,
261            -8.695_404_723_772_031e-4,
262            0.999_999_621_949_602_7,
263            -1.361_752_497_080_270_2e-7,
264            -3.779_734_957_034_089_7e-4,
265            -1.924_880_847_894_457e-7,
266            0.999_999_928_567_997_2,
267        ])
268        .transpose();
269        let act = ReferenceSystem::Iers2010.bias_precession_matrix(time);
270        assert_approx_eq!(act, exp, atol <= 1e-12);
271    }
272
273    #[test]
274    fn test_frame_bias_iau2000() {
275        let exp = DMat3::from_cols_array(&[
276            0.999_999_999_999_994_2,
277            -7.078_279_744_199_197e-8,
278            8.056_217_146_976_134e-8,
279            7.078_279_477_857_338e-8,
280            0.999_999_999_999_997,
281            3.306_041_454_222_136_4e-8,
282            -8.056_217_380_986_972e-8,
283            -3.306_040_883_980_552_3e-8,
284            0.999_999_999_999_996_2,
285        ])
286        .transpose();
287        let act = frame_bias();
288        assert_approx_eq!(act, exp, atol <= 1e-12);
289    }
290
291    #[test]
292    fn test_precession_corrections_iau2000() {
293        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 53736.0);
294        let exp = PrecessionCorrectionsIau2000 {
295            dpsipr: -8.716_465_172_668_348e-8.rad(),
296            depspr: -7.342_018_386_722_813e-9.rad(),
297        };
298        let act = PrecessionCorrectionsIau2000::new(time);
299        assert_approx_eq!(act, exp, atol <= 1e-22);
300    }
301
302    #[test]
303    fn test_decomposition_invariant_iers1996() {
304        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 50123.9999);
305        let bp = ReferenceSystem::Iers1996.bias_precession_matrix(time);
306        let p = ReferenceSystem::Iers1996.precession_matrix(time);
307        let b = frame_bias();
308        assert_approx_eq!(bp, p * b, atol <= 1e-14);
309    }
310
311    #[test]
312    fn test_decomposition_invariant_iers2003() {
313        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 50123.9999);
314        let bp = ReferenceSystem::Iers2003(Iau2000Model::A).bias_precession_matrix(time);
315        let p = ReferenceSystem::Iers2003(Iau2000Model::A).precession_matrix(time);
316        let b = frame_bias();
317        assert_approx_eq!(bp, p * b, atol <= 1e-14);
318    }
319
320    #[test]
321    fn test_decomposition_invariant_iers2010() {
322        let time = Time::from_two_part_julian_date(Tt, 2400000.5, 50123.9999);
323        let bp = ReferenceSystem::Iers2010.bias_precession_matrix(time);
324        let p = ReferenceSystem::Iers2010.precession_matrix(time);
325        let b = frame_bias();
326        assert_approx_eq!(bp, p * b, atol <= 1e-14);
327    }
328}
lox_frames/iers/precession.rs

lox_frames/iers/
precession.rs
