sidereon_core/astro/frames/
precession.rs

1//! IAU 2006 precession matrix and ICRS-to-J2000 frame bias, ported from the
2//! C++ Skyfield-compatible implementation.
3//!
4//! Originally `pub(crate)` inside `orbis_nif`; now public in the core crate so
5//! a Rust-only consumer can reach it without Rustler or the BEAM. The numerics,
6//! summation order, and transcendental sequence are preserved exactly so the
7//! existing Skyfield 0-ULP parity holds.
8//!
9//! All arithmetic uses plain operators (no `f64::mul_add`) so that
10//! rounding matches CPython / Skyfield compiled without FMA contraction.
11
12use crate::astro::constants::time::{DAYS_PER_JULIAN_CENTURY, J2000_JD};
13use crate::astro::constants::units::ARCSEC_TO_RAD;
14use crate::astro::math::mat3::Mat3;
15
16/// Error returned when public precession inputs are outside the valid domain.
17#[derive(Debug, Clone, Copy, PartialEq, Eq, thiserror::Error)]
18pub enum PrecessionError {
19    /// A precession input was non-finite or otherwise invalid.
20    #[error("invalid precession {field}: {reason}")]
21    InvalidInput {
22        field: &'static str,
23        reason: &'static str,
24    },
25}
26
27fn invalid_input(field: &'static str, reason: &'static str) -> PrecessionError {
28    PrecessionError::InvalidInput { field, reason }
29}
30
31fn validate_finite(field: &'static str, value: f64) -> Result<(), PrecessionError> {
32    if value.is_finite() {
33        Ok(())
34    } else {
35        Err(invalid_input(field, "must be finite"))
36    }
37}
38
39fn validate_mat3(field: &'static str, mat: Mat3) -> Result<Mat3, PrecessionError> {
40    if mat.iter().flatten().all(|value| value.is_finite()) {
41        Ok(mat)
42    } else {
43        Err(invalid_input(field, "components must be finite"))
44    }
45}
46
47// ---------------------------------------------------------------------------
48// Precession matrix (IAU 2006)
49// ---------------------------------------------------------------------------
50
51/// Compute the 3x3 precession rotation matrix for the given TDB Julian date,
52/// using the IAU 2006 Fukushima-Williams parameterisation.
53pub fn compute_skyfield_precession_matrix(jd_tdb: f64) -> Result<Mat3, PrecessionError> {
54    validate_finite("jd_tdb", jd_tdb)?;
55    validate_mat3(
56        "precession_matrix",
57        compute_skyfield_precession_matrix_unchecked(jd_tdb),
58    )
59}
60
61pub(crate) fn compute_skyfield_precession_matrix_unchecked(jd_tdb: f64) -> Mat3 {
62    const EPS0_ARCSEC: f64 = 84381.406;
63
64    let t = (jd_tdb - J2000_JD) / DAYS_PER_JULIAN_CENTURY;
65
66    let psia = ((((-0.0000000951 * t + 0.000132851) * t - 0.00114045) * t - 1.0790069) * t
67        + 5038.481507)
68        * t;
69
70    let omegaa =
71        ((((0.0000003337 * t - 0.000000467) * t - 0.00772503) * t + 0.0512623) * t - 0.025754) * t
72            + EPS0_ARCSEC;
73
74    let chia = ((((-0.0000000560 * t + 0.000170663) * t - 0.00121197) * t - 2.3814292) * t
75        + 10.556403)
76        * t;
77
78    let eps0 = EPS0_ARCSEC * ARCSEC_TO_RAD;
79    let psia_rad = psia * ARCSEC_TO_RAD;
80    let omegaa_rad = omegaa * ARCSEC_TO_RAD;
81    let chia_rad = chia * ARCSEC_TO_RAD;
82
83    let sa = eps0.sin();
84    let ca = eps0.cos();
85    let sb = (-psia_rad).sin();
86    let cb = (-psia_rad).cos();
87    let sc = (-omegaa_rad).sin();
88    let cc = (-omegaa_rad).cos();
89    let sd = chia_rad.sin();
90    let cd = chia_rad.cos();
91
92    [
93        [
94            cd * cb - sb * sd * cc,
95            cd * sb * ca + sd * cc * cb * ca - sa * sd * sc,
96            cd * sb * sa + sd * cc * cb * sa + ca * sd * sc,
97        ],
98        [
99            -sd * cb - sb * cd * cc,
100            -sd * sb * ca + cd * cc * cb * ca - sa * cd * sc,
101            -sd * sb * sa + cd * cc * cb * sa + ca * cd * sc,
102        ],
103        [sb * sc, -sc * cb * ca - sa * cc, -sc * cb * sa + cc * ca],
104    ]
105}
106
107// ---------------------------------------------------------------------------
108// ICRS-to-J2000 frame bias matrix
109// ---------------------------------------------------------------------------
110
111/// Build the ICRS-to-J2000 frame bias rotation matrix.
112///
113/// This accounts for the small offset between the ICRS axes and the
114/// mean J2000.0 dynamical frame, parameterised by the three bias angles
115/// xi_0, eta_0, and da_0 from the IAU 2006 conventions.
116pub fn build_icrs_to_j2000() -> Mat3 {
117    let xi0 = -0.0166170 * ARCSEC_TO_RAD;
118    let eta0 = -0.0068192 * ARCSEC_TO_RAD;
119    let da0 = -0.01460 * ARCSEC_TO_RAD;
120
121    let yx = -da0;
122    let zx = xi0;
123    let xy = da0;
124    let zy = eta0;
125    let xz = -xi0;
126    let yz = -eta0;
127
128    [
129        [1.0 - 0.5 * (yx * yx + zx * zx), xy, xz],
130        [yx, 1.0 - 0.5 * (yx * yx + zy * zy), yz],
131        [zx, zy, 1.0 - 0.5 * (zy * zy + zx * zx)],
132    ]
133}
sidereon_core/astro/frames/precession.rs

sidereon_core/astro/frames/
precession.rs
