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

lox_frames/iers/
precession.rs
