|
| 1 | +// SPDX-FileCopyrightText: 2025 Helge Eichhorn <[email protected]> |
| 2 | +// |
| 3 | +// SPDX-License-Identifier: MPL-2.0 |
| 4 | + |
| 5 | +use fast_polynomial::poly_array; |
| 6 | +use lox_test_utils::ApproxEq; |
| 7 | +use lox_time::{Time, julian_dates::JulianDate, time_scales::Tt}; |
| 8 | +use lox_units::Angle; |
| 9 | + |
| 10 | +/// Mean obliquity of the ecliptic. |
| 11 | +#[derive(Debug, Clone, Copy, PartialEq, PartialOrd, ApproxEq)] |
| 12 | +pub struct MeanObliquity(pub Angle); |
| 13 | + |
| 14 | +impl MeanObliquity { |
| 15 | + /// Returns the mean obliquity of the ecliptic based on the IAU1980 precession model. |
| 16 | + pub fn iau1980(time: Time<Tt>) -> Self { |
| 17 | + let t = time.centuries_since_j2000(); |
| 18 | + |
| 19 | + Self(Angle::arcseconds(poly_array( |
| 20 | + t, |
| 21 | + &[84381.448, -46.8150, -0.00059, 0.001813], |
| 22 | + ))) |
| 23 | + } |
| 24 | + |
| 25 | + /// Returns the mean obliquity of the ecliptic based on the IAU2006 precession model. |
| 26 | + pub fn iau2006(time: Time<Tt>) -> Self { |
| 27 | + let t = time.centuries_since_j2000(); |
| 28 | + |
| 29 | + Self(Angle::arcseconds(poly_array( |
| 30 | + t, |
| 31 | + &[ |
| 32 | + 84381.406, |
| 33 | + -46.836769, |
| 34 | + -0.0001831, |
| 35 | + 0.00200340, |
| 36 | + -0.000000576, |
| 37 | + -0.0000000434, |
| 38 | + ], |
| 39 | + ))) |
| 40 | + } |
| 41 | +} |
| 42 | + |
| 43 | +#[cfg(test)] |
| 44 | +mod tests { |
| 45 | + use lox_test_utils::assert_approx_eq; |
| 46 | + |
| 47 | + use super::*; |
| 48 | + |
| 49 | + #[test] |
| 50 | + fn test_mean_obliquity_iau1980() { |
| 51 | + let time = Time::from_two_part_julian_date(Tt, 2400000.5, 54388.0); |
| 52 | + let exp = MeanObliquity(Angle::new(0.409_075_134_764_381_6)); |
| 53 | + let act = MeanObliquity::iau1980(time); |
| 54 | + assert_approx_eq!(act, exp, rtol <= 1e-14); |
| 55 | + } |
| 56 | + |
| 57 | + #[test] |
| 58 | + fn test_mean_obliquity_iau2006() { |
| 59 | + let time = Time::from_two_part_julian_date(Tt, 2400000.5, 54388.0); |
| 60 | + let exp = MeanObliquity(Angle::new(0.409_074_922_938_725_8)); |
| 61 | + let act = MeanObliquity::iau2006(time); |
| 62 | + assert_approx_eq!(act, exp, rtol <= 1e-14); |
| 63 | + } |
| 64 | +} |
0 commit comments