erfa_rust/
G21_safe.rs

1// G21
2//   moon98.c  → eraMoon98_safe
3//   num00a.c  → eraNum00a_safe
4//   num00b.c  → eraNum00b_safe
5//   num06a.c  → eraNum06a_safe
6//   numat.c   → eraNumat_safe
7
8use crate::G19_safe::eraIr_safe;
9use crate::G23_safe::{eraNut06a_safe, eraObl06_safe};
10use crate::G24_safe::eraPfw06_safe;
11use crate::G25_safe::{eraPn00a_safe, eraPn00b_safe};
12use crate::G28_safe::{eraRx_safe, eraRxpv_safe, eraRz_safe};
13use crate::G29_safe::eraS2pv_safe;
14use crate::H1_safe::{ERFA_DAU, ERFA_DD2R, ERFA_DJ00, ERFA_DJC};
15
16#[path = "data/G21_safe/TLR.rs"]
17mod tlr_mod;
18use tlr_mod::TLR;
19
20#[path = "data/G21_safe/TB.rs"]
21mod tb_mod;
22use tb_mod::TB;
23
24pub type ErfaResult<T> = Result<T, ()>;
25
26// eraMoon98_safe: Approximate geocentric Moon pv (Meeus 1998), GCRS, au and au/day.
27pub fn eraMoon98_safe(date1: f64, date2: f64) -> ErfaResult<[[f64; 3]; 2]> {
28    // 1) Fundamental-argument coefficients (degrees).
29    const EL0: [f64; 5] = [
30        218.316_654_36,
31        481_267.881_234_21,
32        -0.001_578_6,
33        1.0 / 538_841.0,
34        -1.0 / 65_194_000.0,
35    ];
36    const D0: [f64; 5] = [
37        297.850_192_1,
38        445_267.111_403_4,
39        -0.001_881_9,
40        1.0 / 545_868.0,
41        1.0 / 113_065_000.0,
42    ];
43    const EM0: [f64; 5] = [
44        357.529_109_2,
45        35_999.050_290_9,
46        -0.000_153_6,
47        1.0 / 24_490_000.0,
48        0.0,
49    ];
50    const EMP0: [f64; 5] = [
51        134.963_396_4,
52        477_198.867_505_5,
53        0.008_741_4,
54        1.0 / 69_699.0,
55        -1.0 / 14_712_000.0,
56    ];
57    const F0: [f64; 5] = [
58        93.272_095_0,
59        483_202.017_523_3,
60        -0.003_653_9,
61        1.0 / 3_526_000.0,
62        1.0 / 863_310_000.0,
63    ];
64
65    // Meeus additional arguments.
66    const A1: (f64, f64) = (119.75, 131.849);
67    const A2: (f64, f64) = (53.09, 479_264.290);
68    const A3: (f64, f64) = (313.45, 481_266.484);
69
70    // Meeus additive-term coefficients.
71    const AL: (f64, f64, f64) = (0.003_958, 0.001_962, 0.000_318);
72    const AB: (f64, f64, f64, f64, f64, f64) = (
73        -0.002_235, 0.000_382, 0.000_175, 0.000_175, 0.000_127, -0.000_115,
74    );
75
76    const R0_METERS: f64 = 385_000_560.0;
77    const E_COEF: (f64, f64) = (-0.002_516, -0.000_007_4);
78
79    // 2) Time argument and fundamental angles.
80    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
81
82    #[inline]
83    fn poly(coeff: &[f64; 5], t: f64) -> (f64, f64) {
84        let w = coeff[0] + (coeff[1] + (coeff[2] + (coeff[3] + coeff[4] * t) * t) * t) * t;
85        let dw = coeff[1] + (2.0 * coeff[2] + (3.0 * coeff[3] + 4.0 * coeff[4] * t) * t) * t;
86        (ERFA_DD2R * (w.rem_euclid(360.0)), ERFA_DD2R * dw)
87    }
88
89    let (elp, delp) = poly(&EL0, t);
90    let (d, dd) = poly(&D0, t);
91    let (em, dem) = poly(&EM0, t);
92    let (emp, demp) = poly(&EMP0, t);
93    let (f, df) = poly(&F0, t);
94
95    let a1 = ERFA_DD2R * (A1.0 + A1.1 * t);
96    let a2 = ERFA_DD2R * (A2.0 + A2.1 * t);
97    let a3 = ERFA_DD2R * (A3.0 + A3.1 * t);
98    let da1 = ERFA_DD2R * AL.0;
99    let da2 = ERFA_DD2R * A2.1;
100    let da3 = ERFA_DD2R * A3.1;
101
102    // E-factor.
103    let e = 1.0 + (E_COEF.0 + E_COEF.1 * t) * t;
104    let de = E_COEF.0 + 2.0 * E_COEF.1 * t;
105    let esq = e * e;
106    let desq = 2.0 * e * de;
107
108    // 3) Meeus additive terms.
109    let elpmf = elp - f;
110    let delpmf = delp - df;
111    let mut vel = AL.0 * a1.sin() + AL.1 * elpmf.sin() + AL.2 * a2.sin();
112    let mut vdel = AL.0 * da1 * a1.cos() + AL.1 * delpmf * elpmf.cos() + AL.2 * da2 * a2.cos();
113
114    let mut vr = 0.0_f64;
115    let mut vdr = 0.0_f64;
116
117    let a1mf = a1 - f;
118    let da1mf = da1 - df;
119    let a1pf = a1 + f;
120    let da1pf = da1 + df;
121    let dlpmp = elp - emp;
122    let slpmp = elp + emp;
123
124    let vb0 = AB.0 * elp.sin()
125        + AB.1 * a3.sin()
126        + AB.2 * a1mf.sin()
127        + AB.3 * a1pf.sin()
128        + AB.4 * dlpmp.sin()
129        + AB.5 * slpmp.sin();
130    let vdb0 = AB.0 * delp * elp.cos()
131        + AB.1 * da3 * a3.cos()
132        + AB.2 * da1mf * a1mf.cos()
133        + AB.3 * da1pf * a1pf.cos()
134        + AB.4 * (delp - demp) * dlpmp.cos()
135        + AB.5 * (delp + demp) * slpmp.cos();
136    let mut vb = vb0;
137    let mut vdb = vdb0;
138
139    // 4) Longitude and distance series.
140    for (nd, nem, nemp, nf, coefl, coefr) in TLR.iter().rev() {
141        let dn = *nd as f64;
142        let emn = *nem as f64;
143        let empn = *nemp as f64;
144        let fn_ = *nf as f64;
145
146        let (en, den) = match nem.abs() {
147            1 => (e, de),
148            2 => (esq, desq),
149            _ => (1.0, 0.0),
150        };
151
152        let arg = dn * d + emn * em + empn * emp + fn_ * f;
153        let darg = dn * dd + emn * dem + empn * demp + fn_ * df;
154
155        let v = arg.sin() * en;
156        let dv = arg.cos() * darg * en + arg.sin() * den;
157        vel += coefl * v;
158        vdel += coefl * dv;
159
160        let v = arg.cos() * en;
161        let dv = -arg.sin() * darg * en + arg.cos() * den;
162        vr += coefr * v;
163        vdr += coefr * dv;
164    }
165    let el = elp + ERFA_DD2R * vel;
166    let del = (delp + ERFA_DD2R * vdel) / ERFA_DJC;
167    let r = (vr + R0_METERS) / ERFA_DAU;
168    let dr = vdr / ERFA_DAU / ERFA_DJC;
169
170    // 5) Latitude series.
171    for (nd, nem, nemp, nf, coefb) in TB.iter().rev() {
172        let dn = *nd as f64;
173        let emn = *nem as f64;
174        let empn = *nemp as f64;
175        let fn_ = *nf as f64;
176        let (en, den) = match nem.abs() {
177            1 => (e, de),
178            2 => (esq, desq),
179            _ => (1.0, 0.0),
180        };
181        let arg = dn * d + emn * em + empn * emp + fn_ * f;
182        let darg = dn * dd + emn * dem + empn * demp + fn_ * df;
183        let v = arg.sin() * en;
184        let dv = arg.cos() * darg * en + arg.sin() * den;
185        vb += coefb * v;
186        vdb += coefb * dv;
187    }
188    let b = vb * ERFA_DD2R;
189    let db = vdb * ERFA_DD2R / ERFA_DJC;
190
191    // 6) Spherical to pv (ecliptic of date), then rotate to GCRS.
192    let pv_local = eraS2pv_safe(el, b, r, del, db, dr)?; // Meeus, ecliptic of date
193
194    // FukushimaWilliams bias-precession, IAU 2006.
195    let (gamb, phib, psib, _epsa) = eraPfw06_safe(date1, date2)?;
196
197    // Build rotation matrix mean ecliptic → GCRS: Rz(psib)  Rx(-phib)  Rz(-gamb).
198    let mut rm = [[0.0_f64; 3]; 3];
199    eraIr_safe(&mut rm)?;
200    eraRz_safe(psib, &mut rm)?;
201    eraRx_safe(-phib, &mut rm)?;
202    eraRz_safe(-gamb, &mut rm)?;
203
204    // Rotate pv into GCRS.
205    let pv_out = eraRxpv_safe(&rm, &pv_local)?;
206    Ok(pv_out)
207}
208
209// eraNum00a_safe: Nutation matrix for IAU 2000A.
210pub fn eraNum00a_safe(date1: f64, date2: f64) -> ErfaResult<[[f64; 3]; 3]> {
211    let (_dpsi, _deps, epsa, _rb, _rp, _rbp, _rn, _rbpn) = eraPn00a_safe(date1, date2)?;
212    let rmatn = eraNumat_safe(epsa, _dpsi, _deps)?;
213    Ok(rmatn)
214}
215
216// eraNum00b_safe: Nutation matrix for IAU 2000B.
217pub fn eraNum00b_safe(date1: f64, date2: f64) -> ErfaResult<[[f64; 3]; 3]> {
218    let (_dpsi, _deps, epsa, _rb, _rp, _rbp, _rn, _rbpn) = eraPn00b_safe(date1, date2)?;
219    let rmatn = eraNumat_safe(epsa, _dpsi, _deps)?;
220    Ok(rmatn)
221}
222
223// eraNum06a_safe: Nutation matrix for IAU 2006/2000A.
224pub fn eraNum06a_safe(date1: f64, date2: f64) -> ErfaResult<[[f64; 3]; 3]> {
225    let eps = eraObl06_safe(date1, date2)?;
226    let (dp, de) = eraNut06a_safe(date1, date2)?;
227    let rmatn = eraNumat_safe(eps, dp, de)?;
228    Ok(rmatn)
229}
230
231// eraNumat_safe: Build nutation matrix from mean obliquity and nutation in longitude/obliquity.
232pub fn eraNumat_safe(epsa: f64, dpsi: f64, deps: f64) -> ErfaResult<[[f64; 3]; 3]> {
233    let mut r = [[0.0_f64; 3]; 3];
234    eraIr_safe(&mut r)?;
235    eraRx_safe(epsa, &mut r)?;
236    eraRz_safe(-dpsi, &mut r)?;
237    eraRx_safe(-(epsa + deps), &mut r)?;
238    Ok(r)
239}
erfa_rust/G21_safe.rs

erfa_rust/
G21_safe.rs
