erfa_rust/
G22_safe.rs

1// G22
2//   nut00a.c → eraNut00a_safe
3
4use crate::G15_safe::{
5    eraFae03_safe, eraFaf03_safe, eraFaju03_safe, eraFal03_safe, eraFama03_safe, eraFame03_safe,
6    eraFaom03_safe, eraFapa03_safe, eraFasa03_safe, eraFaur03_safe, eraFave03_safe,
7};
8use crate::H1_safe::{ERFA_D2PI, ERFA_DAS2R, ERFA_DJ00, ERFA_DJC, ERFA_TURNAS};
9
10#[path = "data/G22_safe/XLS.rs"]
11mod xls_mod;
12use xls_mod::XLS;
13
14#[path = "data/G22_safe/XPL.rs"]
15mod xpl_mod;
16use xpl_mod::XPL;
17
18pub type ErfaResult<T> = Result<T, ()>;
19
20// Small helper: positive modulo for floating angles.
21#[inline]
22fn fmod_pos(a: f64, p: f64) -> f64 {
23    (a % p + p) % p
24}
25
26// coefficient tables
27// Types must match the data tables defined in data/G22/*.rs (verbatim).
28#[derive(Clone, Copy)]
29struct LuniSolar {
30    nl: i32,
31    nlp: i32,
32    nf: i32,
33    nd: i32,
34    nom: i32,
35    sp: f64,
36    spt: f64,
37    cp: f64,
38    ce: f64,
39    cet: f64,
40    se: f64,
41}
42
43#[derive(Clone, Copy)]
44struct Planetary {
45    nl: i32,
46    nf: i32,
47    nd: i32,
48    nom: i32,
49    nme: i32,
50    nve: i32,
51    nea: i32,
52    nma: i32,
53    nju: i32,
54    nsa: i32,
55    nur: i32,
56    nne: i32,
57    npa: i32,
58    sp: i32,
59    cp: i32,
60    se: i32,
61    ce: i32,
62}
63
64// main routine
65// Nutation, IAU-2000A model (no free-core nutation)
66pub fn eraNut00a_safe(date1: f64, date2: f64) -> ErfaResult<(f64, f64)> {
67    // 0.1 µas to radians.
68    const U2R: f64 = ERFA_DAS2R / 1e7;
69
70    // Time in Julian centuries from J2000.0.
71    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
72
73    // Luni-solar fundamental arguments.
74    let el = eraFal03_safe(t)?; // Moon mean anomaly
75    // Sun mean anomaly (MHB2000 explicit series in arcsec, then radians).
76    let elp = fmod_pos(
77        1_287_104.793_05
78            + t * (129_596_581.0481 + t * (-0.5532 + t * (0.000_136 + t * (-0.000_011_49)))),
79        ERFA_TURNAS,
80    ) * ERFA_DAS2R;
81    let f = eraFaf03_safe(t)?; // Moon argument of latitude
82    let d = fmod_pos(
83        1_072_260.703_69
84            + t * (1_602_961_601.2090 + t * (-6.3706 + t * (0.006_593 + t * (-0.000_031_69)))),
85        ERFA_TURNAS,
86    ) * ERFA_DAS2R;
87    let om = eraFaom03_safe(t)?; // Moon ascending node longitude
88
89    // Accumulators for luni-solar series.
90    let mut dp = 0.0f64;
91    let mut de = 0.0f64;
92
93    // Sum luni-solar terms (reverse order to match C loop order).
94    for term in XLS.iter().rev() {
95        let arg = fmod_pos(
96            term.nl as f64 * el
97                + term.nlp as f64 * elp
98                + term.nf as f64 * f
99                + term.nd as f64 * d
100                + term.nom as f64 * om,
101            ERFA_D2PI,
102        );
103        let s = arg.sin();
104        let c = arg.cos();
105
106        dp += (term.sp + term.spt * t) * s + term.cp * c;
107        de += (term.ce + term.cet * t) * c + term.se * s;
108    }
109    let dpsils = dp * U2R;
110    let depsls = de * U2R;
111
112    // Planetary contributions (MHB2000).
113    let al = fmod_pos(2.355_555_98 + 8_328.691_426_9554 * t, ERFA_D2PI);
114    let af = fmod_pos(1.627_905_234 + 8_433.466_158_1310 * t, ERFA_D2PI);
115    let ad = fmod_pos(5.198_466_741 + 7_771.377_146_8121 * t, ERFA_D2PI);
116    let aom = fmod_pos(2.182_439_20 - 33.757_045 * t, ERFA_D2PI);
117    let apa = eraFapa03_safe(t)?; // General precession in longitude
118
119    // Planetary longitudes.
120    let alme = eraFame03_safe(t)?;
121    let alve = eraFave03_safe(t)?;
122    let alea = eraFae03_safe(t)?;
123    let alma = eraFama03_safe(t)?;
124    let alju = eraFaju03_safe(t)?;
125    let alsa = eraFasa03_safe(t)?;
126    let alur = eraFaur03_safe(t)?;
127    let alne = fmod_pos(5.321_159_000 + 3.812_777_4000 * t, ERFA_D2PI);
128
129    // Reset accumulators for planetary series.
130    dp = 0.0;
131    de = 0.0;
132
133    for term in XPL.iter().rev() {
134        let arg = fmod_pos(
135            term.nl as f64 * al
136                + term.nf as f64 * af
137                + term.nd as f64 * ad
138                + term.nom as f64 * aom
139                + term.nme as f64 * alme
140                + term.nve as f64 * alve
141                + term.nea as f64 * alea
142                + term.nma as f64 * alma
143                + term.nju as f64 * alju
144                + term.nsa as f64 * alsa
145                + term.nur as f64 * alur
146                + term.nne as f64 * alne
147                + term.npa as f64 * apa,
148            ERFA_D2PI,
149        );
150        let s = arg.sin();
151        let c = arg.cos();
152
153        dp += term.sp as f64 * s + term.cp as f64 * c;
154        de += term.se as f64 * s + term.ce as f64 * c;
155    }
156    let dpsipl = dp * U2R;
157    let depspl = de * U2R;
158
159    // Return total nutation in longitude and obliquity.
160    Ok((dpsils + dpsipl, depsls + depspl))
161}
erfa_rust/G22_safe.rs

erfa_rust/
G22_safe.rs
