erfa_rust/
G24_safe.rs

1// G24
2//   p06e.c   → eraP06e_safe
3//   p2pv.c   → eraP2pv_safe
4//   p2s.c    → eraP2s_safe
5//   pap.c    → eraPap_safe
6//   pas.c    → eraPas_safe
7//   pb06.c   → eraPb06_safe
8//   pdp.c    → eraPdp_safe
9//   pfw06.c  → eraPfw06_safe
10//   plan94.c → eraPlan94_safe
11//   pm.c     → eraPm_safe
12
13use crate::G1_safe::eraAnpm_safe;
14use crate::G23_safe::eraObl06_safe;
15use crate::G25_safe::{eraPmat06_safe, eraPmp_safe, eraPn_safe};
16use crate::G27_safe::eraPxp_safe;
17use crate::G28_safe::eraRz_safe;
18use crate::G7_safe::eraC2s_safe;
19use crate::H1_safe::{ERFA_D2PI, ERFA_DAS2R, ERFA_DJ00, ERFA_DJC, ERFA_DJM};
20
21#[path = "data/G24_safe/plan94_tables.rs"]
22mod plan94_tables;
23use plan94_tables::{
24    ASCENDING_NODE as OMEGA, ECCENTRICITY as E, INCLINATION as DINC, MEAN_LONGITUDE as DLM,
25    PERIHELION as PI_, PLANET_MASSES as AMAS, SEMI_MAJOR_AXIS as A, TRIG_ARG_LONGITUDE as KQ,
26    TRIG_ARG_PERIHELION as KP, TRIG_COEFF_A_COS as CA, TRIG_COEFF_A_SIN as SA,
27    TRIG_COEFF_L_COS as CL, TRIG_COEFF_L_SIN as SL,
28};
29
30pub type ErfaResult<T> = Result<T, ()>;
31
32// eraP06e_safe: IAU 2006 equinox-based precession parameter set.
33pub fn eraP06e_safe(
34    date1: f64,
35    date2: f64,
36) -> ErfaResult<(
37    f64,
38    f64,
39    f64,
40    f64,
41    f64,
42    f64,
43    f64,
44    f64,
45    f64,
46    f64,
47    f64,
48    f64,
49    f64,
50    f64,
51    f64,
52    f64,
53)> {
54    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
55
56    let eps0 = 84381.406 * ERFA_DAS2R;
57
58    let psia = (5038.481_507
59        + (-1.079_006_9 + (-0.001_140_45 + (0.000_132_851 + (-0.000_000_0951) * t) * t) * t) * t)
60        * t
61        * ERFA_DAS2R;
62
63    let oma = eps0
64        + ((-0.025_754
65            + (0.051_262_3 + (-0.007_725_03 + (-0.000_000_467 + 0.000_000_3337 * t) * t) * t) * t)
66            * t)
67            * ERFA_DAS2R;
68
69    let bpa = (4.199_094
70        + (0.193_987_3 + (-0.000_224_66 + (-0.000_000_912 + 0.000_000_0120 * t) * t) * t) * t)
71        * t
72        * ERFA_DAS2R;
73
74    let bqa = (-46.811_015
75        + (0.051_028_3 + (0.000_524_13 + (-0.000_000_646 + (-0.000_000_0172) * t) * t) * t) * t)
76        * t
77        * ERFA_DAS2R;
78
79    let pia = (46.998_973
80        + (-0.033_492_6 + (-0.000_125_59 + (0.000_000_113 + (-0.000_000_0022) * t) * t) * t) * t)
81        * t
82        * ERFA_DAS2R;
83
84    let bpia = (629_546.7936
85        + (-867.95758
86            + (0.157_992 + (-0.000_5371 + (-0.000_047_97 + 0.000_000_072 * t) * t) * t) * t)
87            * t)
88        * ERFA_DAS2R;
89
90    let epsa = eraObl06_safe(date1, date2)?;
91
92    let chia = (10.556_403
93        + (-2.381_429_2 + (-0.001_211_97 + (0.000_170_663 + (-0.000_000_0560) * t) * t) * t) * t)
94        * t
95        * ERFA_DAS2R;
96
97    let za = (-2.650_545
98        + (2306.077_181
99            + (1.092_734_8 + (0.018_268_37 + (-0.000_028_596 + (-0.000_000_2904) * t) * t) * t)
100                * t)
101            * t)
102        * ERFA_DAS2R;
103
104    let zetaa = (2.650_545
105        + (2306.083_227
106            + (0.298_849_9 + (0.018_018_28 + (-0.000_005_971 + (-0.000_000_3173) * t) * t) * t)
107                * t)
108            * t)
109        * ERFA_DAS2R;
110
111    let thetaa = (2004.191_903
112        + (-0.429_493_4 + (-0.041_822_64 + (-0.000_007_089 + (-0.000_000_1274) * t) * t) * t) * t)
113        * t
114        * ERFA_DAS2R;
115
116    let pa = (5028.796_195
117        + (1.105_434_8 + (0.000_079_64 + (-0.000_023_857 + (-0.000_000_0383) * t) * t) * t) * t)
118        * t
119        * ERFA_DAS2R;
120
121    let gam = (10.556_403
122        + (0.493_204_4 + (-0.000_312_38 + (-0.000_002_788 + 0.000_000_0260 * t) * t) * t) * t)
123        * t
124        * ERFA_DAS2R;
125
126    let phi = eps0
127        + ((-46.811_015
128            + (0.051_126_9 + (0.000_532_89 + (-0.000_000_440 + (-0.000_000_0176) * t) * t) * t)
129                * t)
130            * t)
131            * ERFA_DAS2R;
132
133    let psi = (5038.481_507
134        + (1.558_417_6 + (-0.000_185_22 + (-0.000_026_452 + (-0.000_000_0148) * t) * t) * t) * t)
135        * t
136        * ERFA_DAS2R;
137
138    Ok((
139        eps0, psia, oma, bpa, bqa, pia, bpia, epsa, chia, za, zetaa, thetaa, pa, gam, phi, psi,
140    ))
141}
142
143// eraP2pv_safe: Extend p-vector to pv-vector with zero velocity.
144pub fn eraP2pv_safe(p: &[f64; 3]) -> ErfaResult<[[f64; 3]; 2]> {
145    let mut pv = [[0.0_f64; 3]; 2];
146    pv[0] = *p;
147    Ok(pv)
148}
149
150// eraP2s_safe: Convert p-vector to spherical coordinates (theta, phi, r).
151pub fn eraP2s_safe(p: &[f64; 3]) -> ErfaResult<(f64, f64, f64)> {
152    let (theta, phi) = eraC2s_safe(p)?;
153    let r = eraPm_safe(p)?;
154    Ok((theta, phi, r))
155}
156
157// eraPap_safe: Position angle from two p-vectors (radians).
158pub fn eraPap_safe(a: &[f64; 3], b: &[f64; 3]) -> ErfaResult<f64> {
159    let (am, au) = eraPn_safe(a)?;
160    let bm = eraPm_safe(b)?;
161    let (st, ct) = if am == 0.0 || bm == 0.0 {
162        (0.0, 1.0)
163    } else {
164        let xa = a[0];
165        let ya = a[1];
166        let za = a[2];
167        let eta = [-xa * za, -ya * za, xa * xa + ya * ya];
168        let xi = eraPxp_safe(&eta, &au)?;
169        let a2b = eraPmp_safe(b, a)?;
170        let st_val = eraPdp_safe(&a2b, &xi)?;
171        let mut ct_val = eraPdp_safe(&a2b, &eta)?;
172        if st_val == 0.0 && ct_val == 0.0 {
173            ct_val = 1.0;
174        }
175        (st_val, ct_val)
176    };
177    Ok(st.atan2(ct))
178}
179
180// eraPas_safe: Position angle from two spherical positions.
181pub fn eraPas_safe(al: f64, ap: f64, bl: f64, bp: f64) -> ErfaResult<f64> {
182    let dl = bl - al;
183    let y = dl.sin() * bp.cos();
184    let x = bp.sin() * ap.cos() - bp.cos() * ap.sin() * dl.cos();
185    Ok(if x != 0.0 || y != 0.0 {
186        y.atan2(x)
187    } else {
188        0.0
189    })
190}
191
192// eraPb06_safe: Precession-bias Euler angles (bzeta, bz, btheta), IAU 2006.
193pub fn eraPb06_safe(date1: f64, date2: f64) -> ErfaResult<(f64, f64, f64)> {
194    let r = eraPmat06_safe(date1, date2)?;
195
196    let mut y = r[1][2];
197    let mut x = -r[0][2];
198    let bz = if x < 0.0 {
199        y = -y;
200        x = -x;
201        if x != 0.0 || y != 0.0 {
202            -y.atan2(x)
203        } else {
204            0.0
205        }
206    } else if x != 0.0 || y != 0.0 {
207        -y.atan2(x)
208    } else {
209        0.0
210    };
211
212    let mut r2 = r;
213    eraRz_safe(bz, &mut r2)?;
214
215    y = r2[0][2];
216    x = r2[2][2];
217    let btheta = if x != 0.0 || y != 0.0 {
218        -y.atan2(x)
219    } else {
220        0.0
221    };
222
223    y = -r2[1][0];
224    x = r2[1][1];
225    let bzeta = if x != 0.0 || y != 0.0 {
226        -y.atan2(x)
227    } else {
228        0.0
229    };
230
231    Ok((bzeta, bz, btheta))
232}
233
234// eraPdp_safe: Dot product of two p-vectors.
235pub fn eraPdp_safe(a: &[f64; 3], b: &[f64; 3]) -> ErfaResult<f64> {
236    Ok(a[0] * b[0] + a[1] * b[1] + a[2] * b[2])
237}
238
239// eraPfw06_safe: FukushimaWilliams angles (gamb, phib, psib, epsa), IAU 2006.
240pub fn eraPfw06_safe(date1: f64, date2: f64) -> ErfaResult<(f64, f64, f64, f64)> {
241    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
242
243    let gamb = (-0.052_928
244        + (10.556_378
245            + (0.493_204_4 + (-0.000_312_38 + (-0.000_002_788 + 0.000_000_0260 * t) * t) * t) * t)
246            * t)
247        * ERFA_DAS2R;
248
249    let phib = (84_381.412_819
250        + (-46.811_016
251            + (0.051_126_8 + (0.000_532_89 + (-0.000_000_440 + (-0.000_000_0176) * t) * t) * t)
252                * t)
253            * t)
254        * ERFA_DAS2R;
255
256    let psib = (-0.041_775
257        + (5038.481_484
258            + (1.558_417_5 + (-0.000_185_22 + (-0.000_026_452 + (-0.000_000_0148) * t) * t) * t)
259                * t)
260            * t)
261        * ERFA_DAS2R;
262
263    let epsa = eraObl06_safe(date1, date2)?;
264
265    Ok((gamb, phib, psib, epsa))
266}
267
268// eraPm_safe: Modulus of a 3-vector.
269pub fn eraPm_safe(p: &[f64; 3]) -> ErfaResult<f64> {
270    Ok((p[0] * p[0] + p[1] * p[1] + p[2] * p[2]).sqrt())
271}
272
273// eraPlan94_safe: Approximate heliocentric position/velocity for planet np (1..8).
274pub fn eraPlan94_safe(date1: f64, date2: f64, np_in: i32) -> ErfaResult<([[f64; 3]; 2], i32)> {
275    const GK: f64 = 0.017_202_098_950; // Gaussian gravitational constant (au^3/d^2).
276    const SINEPS: f64 = 0.397_777_155_931_913_7; // J2000 mean obliquity (IAU 1976).
277    const COSEPS: f64 = 0.917_482_062_069_181_8;
278    const KMAX: i32 = 10;
279
280    if np_in < 1 || np_in > 8 {
281        return Ok(([[0.0; 3]; 2], -1));
282    }
283    let np = (np_in - 1) as usize;
284
285    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJM;
286    let mut jstat = if t.abs() <= 1.0 { 0 } else { 1 };
287
288    let mut da = A[np][0] + (A[np][1] + A[np][2] * t) * t;
289    let mut dl = (3600.0 * DLM[np][0] + (DLM[np][1] + DLM[np][2] * t) * t) * ERFA_DAS2R;
290    let de = E[np][0] + (E[np][1] + E[np][2] * t) * t;
291    let dp = eraAnpm_safe((3600.0 * PI_[np][0] + (PI_[np][1] + PI_[np][2] * t) * t) * ERFA_DAS2R)?;
292    let di = (3600.0 * DINC[np][0] + (DINC[np][1] + DINC[np][2] * t) * t) * ERFA_DAS2R;
293    let dom =
294        eraAnpm_safe((3600.0 * OMEGA[np][0] + (OMEGA[np][1] + OMEGA[np][2] * t) * t) * ERFA_DAS2R)?;
295
296    let dmu = 0.359_536_20 * t;
297    for k in 0..8 {
298        let arga = KP[np][k] as f64 * dmu;
299        let argl = KQ[np][k] as f64 * dmu;
300        da += (CA[np][k] as f64 * arga.cos() + SA[np][k] as f64 * arga.sin()) * 1e-7;
301        dl += (CL[np][k] as f64 * argl.cos() + SL[np][k] as f64 * argl.sin()) * 1e-7;
302    }
303    let arga = KP[np][8] as f64 * dmu;
304    da += t * (CA[np][8] as f64 * arga.cos() + SA[np][8] as f64 * arga.sin()) * 1e-7;
305    for k in 8..10 {
306        let argl = KQ[np][k] as f64 * dmu;
307        dl += t * (CL[np][k] as f64 * argl.cos() + SL[np][k] as f64 * argl.sin()) * 1e-7;
308    }
309    dl = dl.rem_euclid(ERFA_D2PI);
310
311    let am = dl - dp;
312    let mut ae = am + de * am.sin();
313    let mut k_iter = 0;
314    loop {
315        let dae = (am - ae + de * ae.sin()) / (1.0 - de * ae.cos());
316        ae += dae;
317        k_iter += 1;
318        if k_iter >= KMAX || dae.abs() <= 1e-12 {
319            break;
320        }
321    }
322    if k_iter >= KMAX {
323        jstat = 2;
324    }
325
326    let ae2 = 0.5 * ae;
327    let at = 2.0 * ((((1.0 + de) / (1.0 - de)).sqrt()) * ae2.sin()).atan2(ae2.cos());
328
329    let r = da * (1.0 - de * ae.cos());
330    let v = GK * ((1.0 + 1.0 / AMAS[np]) / (da * da * da)).sqrt();
331
332    let si2 = (0.5 * di).sin();
333    let xq = si2 * dom.cos();
334    let xp = si2 * dom.sin();
335    let tl = at + dp;
336    let (xsw, xcw) = tl.sin_cos();
337    let xm2 = 2.0 * (xp * xcw - xq * xsw);
338    let xf = da / (1.0 - de * de).sqrt();
339    let ci2 = (0.5 * di).cos();
340    let xms = (de * dp.sin() + xsw) * xf;
341    let xmc = (de * dp.cos() + xcw) * xf;
342    let xpxq2 = 2.0 * xp * xq;
343
344    let x = r * (xcw - xm2 * xp);
345    let y = r * (xsw + xm2 * xq);
346    let z = r * (-xm2 * ci2);
347
348    let mut pv = [[0.0_f64; 3]; 2];
349    pv[0][0] = x;
350    pv[0][1] = y * COSEPS - z * SINEPS;
351    pv[0][2] = y * SINEPS + z * COSEPS;
352
353    let xdot = v * ((-1.0 + 2.0 * xp * xp) * xms + xpxq2 * xmc);
354    let ydot = v * ((1.0 - 2.0 * xq * xq) * xmc - xpxq2 * xms);
355    let zdot = v * (2.0 * ci2 * (xp * xms + xq * xmc));
356
357    pv[1][0] = xdot;
358    pv[1][1] = ydot * COSEPS - zdot * SINEPS;
359    pv[1][2] = ydot * SINEPS + zdot * COSEPS;
360
361    Ok((pv, jstat))
362}
erfa_rust/G24_safe.rs

erfa_rust/
G24_safe.rs
