erfa_rust/
G11_safe.rs

1// G11
2//   eceq06.c  → eraEceq06_safe
3//   ecm06.c   → eraEcm06_safe
4//   ee00.c    → eraEe00_safe
5//   ee00a.c   → eraEe00a_safe
6//   ee00b.c   → eraEe00b_safe
7//   ee06a.c   → eraEe06a_safe
8//   eect00.c  → eraEect00_safe
9//   eform.c   → eraEform_safe
10//   eo06a.c   → eraEo06a_safe
11//   eors.c    → eraEors_safe
12
13use crate::G15_safe::{
14    eraFad03_safe, eraFae03_safe, eraFaf03_safe, eraFal03_safe, eraFalp03_safe, eraFaom03_safe,
15    eraFapa03_safe, eraFave03_safe,
16};
17use crate::G17_safe::{eraGmst06_safe, eraGst06a_safe};
18use crate::G19_safe::eraIr_safe;
19use crate::G1_safe::{eraAnp_safe, eraAnpm_safe};
20use crate::G22_safe::eraNut00a_safe;
21use crate::G23_safe::{eraNut00b_safe, eraObl06_safe, eraObl80_safe};
22use crate::G25_safe::eraPmat06_safe;
23use crate::G26_safe::{eraPnm06a_safe, eraPr00_safe};
24use crate::G28_safe::{eraRx_safe, eraRxr_safe};
25use crate::G29_safe::eraS2c_safe;
26use crate::G30_safe::eraS06_safe;
27use crate::G33_safe::eraTrxp_safe;
28use crate::G6_safe::eraBpn2xy_safe;
29use crate::G7_safe::eraC2s_safe;
30use crate::H1_safe::{ERFA_DAS2R, ERFA_DJ00, ERFA_DJC, ERFA_GRS80, ERFA_WGS72, ERFA_WGS84};
31
32pub type ErfaResult<T> = Result<T, ()>;
33
34// Ecliptic to ICRS (equatorial), IAU 2006.
35pub fn eraEceq06_safe(date1: f64, date2: f64, dl: f64, db: f64) -> ErfaResult<(f64, f64)> {
36    let v1 = eraS2c_safe(dl, db)?;
37    let mut rm = [[0.0_f64; 3]; 3];
38    eraEcm06_safe(date1, date2, &mut rm)?;
39    let v2 = eraTrxp_safe(&rm, &v1)?;
40    let (a, b) = eraC2s_safe(&v2)?;
41    let dr = eraAnp_safe(a)?;
42    let dd = eraAnpm_safe(b)?;
43    Ok((dr, dd))
44}
45
46// Rotation matrix, ICRS equatorial to ecliptic (IAU 2006).
47pub fn eraEcm06_safe(date1: f64, date2: f64, rm: &mut [[f64; 3]; 3]) -> ErfaResult<()> {
48    let ob = eraObl06_safe(date1, date2)?;
49    let bp = eraPmat06_safe(date1, date2)?;
50    let mut e = [[0.0_f64; 3]; 3];
51    eraIr_safe(&mut e)?;
52    eraRx_safe(ob, &mut e)?;
53    let m = eraRxr_safe(&e, &bp)?;
54    *rm = m;
55    Ok(())
56}
57
58// Equation of the equinoxes (IAU 2000), given epsa and dpsi plus complementary terms.
59pub fn eraEe00_safe(date1: f64, date2: f64, epsa: f64, dpsi: f64) -> ErfaResult<f64> {
60    let ee = dpsi * epsa.cos() + eraEect00_safe(date1, date2)?;
61    Ok(ee)
62}
63
64// Equation of the equinoxes consistent with IAU 2000A nutation.
65pub fn eraEe00a_safe(date1: f64, date2: f64) -> ErfaResult<f64> {
66    let (_dpsipr, depspr) = eraPr00_safe(date1, date2)?;
67    let epsa = eraObl80_safe(date1, date2)? + depspr;
68    let (dpsi, _deps) = eraNut00a_safe(date1, date2)?;
69    eraEe00_safe(date1, date2, epsa, dpsi)
70}
71
72// Equation of the equinoxes consistent with IAU 2000B nutation.
73pub fn eraEe00b_safe(date1: f64, date2: f64) -> ErfaResult<f64> {
74    let (_dpsipr, depspr) = eraPr00_safe(date1, date2)?;
75    let epsa = eraObl80_safe(date1, date2)? + depspr;
76    let (dpsi, _deps) = eraNut00b_safe(date1, date2)?;
77    eraEe00_safe(date1, date2, epsa, dpsi)
78}
79
80// Equation of the equinoxes, IAU 2006/2000A.
81pub fn eraEe06a_safe(date1: f64, date2: f64) -> ErfaResult<f64> {
82    let gst06a = eraGst06a_safe(0.0, 0.0, date1, date2)?;
83    let gmst06 = eraGmst06_safe(0.0, 0.0, date1, date2)?;
84    let ee = eraAnpm_safe(gst06a - gmst06)?;
85    Ok(ee)
86}
87
88// Complementary terms for the equation of the equinoxes (IAU 2000), radians.
89pub fn eraEect00_safe(date1: f64, date2: f64) -> ErfaResult<f64> {
90    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
91
92    #[repr(C)]
93    struct TERM {
94        nfa: [i32; 8],
95        s: f64,
96        c: f64,
97    }
98
99    // Terms of order t^0
100    const E0: &[TERM] = &[
101        TERM {
102            nfa: [0, 0, 0, 0, 1, 0, 0, 0],
103            s: 2640.96e-6,
104            c: -0.39e-6,
105        },
106        TERM {
107            nfa: [0, 0, 0, 0, 2, 0, 0, 0],
108            s: 63.52e-6,
109            c: -0.02e-6,
110        },
111        TERM {
112            nfa: [0, 0, 2, -2, 3, 0, 0, 0],
113            s: 11.75e-6,
114            c: 0.01e-6,
115        },
116        TERM {
117            nfa: [0, 0, 2, -2, 1, 0, 0, 0],
118            s: 11.21e-6,
119            c: 0.01e-6,
120        },
121        TERM {
122            nfa: [0, 0, 2, -2, 2, 0, 0, 0],
123            s: -4.55e-6,
124            c: 0.00e-6,
125        },
126        TERM {
127            nfa: [0, 0, 2, 0, 3, 0, 0, 0],
128            s: 2.02e-6,
129            c: 0.00e-6,
130        },
131        TERM {
132            nfa: [0, 0, 2, 0, 1, 0, 0, 0],
133            s: 1.98e-6,
134            c: 0.00e-6,
135        },
136        TERM {
137            nfa: [0, 0, 0, 0, 3, 0, 0, 0],
138            s: -1.72e-6,
139            c: 0.00e-6,
140        },
141        TERM {
142            nfa: [0, 1, 0, 0, 1, 0, 0, 0],
143            s: -1.41e-6,
144            c: -0.01e-6,
145        },
146        TERM {
147            nfa: [0, 1, 0, 0, -1, 0, 0, 0],
148            s: -1.26e-6,
149            c: -0.01e-6,
150        },
151        TERM {
152            nfa: [1, 0, 0, 0, -1, 0, 0, 0],
153            s: -0.63e-6,
154            c: 0.00e-6,
155        },
156        TERM {
157            nfa: [1, 0, 0, 0, 1, 0, 0, 0],
158            s: -0.63e-6,
159            c: 0.00e-6,
160        },
161        TERM {
162            nfa: [0, 1, 2, -2, 3, 0, 0, 0],
163            s: 0.46e-6,
164            c: 0.00e-6,
165        },
166        TERM {
167            nfa: [0, 1, 2, -2, 1, 0, 0, 0],
168            s: 0.45e-6,
169            c: 0.00e-6,
170        },
171        TERM {
172            nfa: [0, 0, 4, -4, 4, 0, 0, 0],
173            s: 0.36e-6,
174            c: 0.00e-6,
175        },
176        TERM {
177            nfa: [0, 0, 1, -1, 1, -8, 12, 0],
178            s: -0.24e-6,
179            c: -0.12e-6,
180        },
181        TERM {
182            nfa: [0, 0, 2, 0, 0, 0, 0, 0],
183            s: 0.32e-6,
184            c: 0.00e-6,
185        },
186        TERM {
187            nfa: [0, 0, 2, 0, 2, 0, 0, 0],
188            s: 0.28e-6,
189            c: 0.00e-6,
190        },
191        TERM {
192            nfa: [1, 0, 2, 0, 3, 0, 0, 0],
193            s: 0.27e-6,
194            c: 0.00e-6,
195        },
196        TERM {
197            nfa: [1, 0, 2, 0, 1, 0, 0, 0],
198            s: 0.26e-6,
199            c: 0.00e-6,
200        },
201        TERM {
202            nfa: [0, 0, 2, -2, 0, 0, 0, 0],
203            s: -0.21e-6,
204            c: 0.00e-6,
205        },
206        TERM {
207            nfa: [0, 1, -2, 2, -3, 0, 0, 0],
208            s: 0.19e-6,
209            c: 0.00e-6,
210        },
211        TERM {
212            nfa: [0, 1, -2, 2, -1, 0, 0, 0],
213            s: 0.18e-6,
214            c: 0.00e-6,
215        },
216        TERM {
217            nfa: [0, 0, 0, 0, 0, 8, -13, -1],
218            s: -0.10e-6,
219            c: 0.05e-6,
220        },
221        TERM {
222            nfa: [0, 0, 0, 2, 0, 0, 0, 0],
223            s: 0.15e-6,
224            c: 0.00e-6,
225        },
226        TERM {
227            nfa: [2, 0, -2, 0, -1, 0, 0, 0],
228            s: -0.14e-6,
229            c: 0.00e-6,
230        },
231        TERM {
232            nfa: [1, 0, 0, -2, 1, 0, 0, 0],
233            s: 0.14e-6,
234            c: 0.00e-6,
235        },
236        TERM {
237            nfa: [0, 1, 2, -2, 2, 0, 0, 0],
238            s: -0.14e-6,
239            c: 0.00e-6,
240        },
241        TERM {
242            nfa: [1, 0, 0, -2, -1, 0, 0, 0],
243            s: 0.14e-6,
244            c: 0.00e-6,
245        },
246        TERM {
247            nfa: [0, 0, 4, -2, 4, 0, 0, 0],
248            s: 0.13e-6,
249            c: 0.00e-6,
250        },
251        TERM {
252            nfa: [0, 0, 2, -2, 4, 0, 0, 0],
253            s: -0.11e-6,
254            c: 0.00e-6,
255        },
256        TERM {
257            nfa: [1, 0, -2, 0, -3, 0, 0, 0],
258            s: 0.11e-6,
259            c: 0.00e-6,
260        },
261        TERM {
262            nfa: [1, 0, -2, 0, -1, 0, 0, 0],
263            s: 0.11e-6,
264            c: 0.00e-6,
265        },
266    ];
267
268    // Terms of order t^1
269    const E1: &[TERM] = &[TERM {
270        nfa: [0, 0, 0, 0, 1, 0, 0, 0],
271        s: -0.87e-6,
272        c: 0.00e-6,
273    }];
274
275    const NE0: usize = 33;
276    const NE1: usize = 1;
277
278    let mut fa = [0.0_f64; 14];
279    fa[0] = eraFal03_safe(t)?;
280    fa[1] = eraFalp03_safe(t)?;
281    fa[2] = eraFaf03_safe(t)?;
282    fa[3] = eraFad03_safe(t)?;
283    fa[4] = eraFaom03_safe(t)?;
284    fa[5] = eraFave03_safe(t)?;
285    fa[6] = eraFae03_safe(t)?;
286    fa[7] = eraFapa03_safe(t)?;
287
288    let mut s0 = 0.0_f64;
289    let mut s1 = 0.0_f64;
290
291    for i in (0..NE0).rev() {
292        let mut a = 0.0_f64;
293        for j in 0..8 {
294            a += (E0[i].nfa[j] as f64) * fa[j];
295        }
296        s0 += E0[i].s * a.sin() + E0[i].c * a.cos();
297    }
298
299    for i in (0..NE1).rev() {
300        let mut a = 0.0_f64;
301        for j in 0..8 {
302            a += (E1[i].nfa[j] as f64) * fa[j];
303        }
304        s1 += E1[i].s * a.sin() + E1[i].c * a.cos();
305    }
306
307    let eect = (s0 + s1 * t) * ERFA_DAS2R;
308    Ok(eect)
309}
310
311// Reference ellipsoid parameters; returns ((a, f), j).
312pub fn eraEform_safe(n: i32) -> ErfaResult<((f64, f64), i32)> {
313    let (a, f, j) = match n {
314        ERFA_WGS84 => (6_378_137.0, 1.0 / 298.257_223_563, 0),
315        ERFA_GRS80 => (6_378_137.0, 1.0 / 298.257_222_101, 0),
316        ERFA_WGS72 => (6_378_135.0, 1.0 / 298.26, 0),
317        _ => (0.0, 0.0, -1),
318    };
319    Ok(((a, f), j))
320}
321
322// Equation of origins, IAU 2006/2000A.
323pub fn eraEo06a_safe(date1: f64, date2: f64) -> ErfaResult<f64> {
324    let r = eraPnm06a_safe(date1, date2)?;
325    let (x, y) = eraBpn2xy_safe(&r)?;
326    let s = eraS06_safe(date1, date2, x, y)?;
327    let eo = eraEors_safe(&r, s)?;
328    Ok(eo)
329}
330
331// Equation of the origins, given NPB matrix and s.
332pub fn eraEors_safe(rnpb: &[[f64; 3]; 3], s: f64) -> ErfaResult<f64> {
333    let x = rnpb[2][0];
334    let ax = x / (1.0 + rnpb[2][2]);
335    let xs = 1.0 - ax * x;
336    let ys = -ax * rnpb[2][1];
337    let zs = -x;
338    let p = rnpb[0][0] * xs + rnpb[0][1] * ys + rnpb[0][2] * zs;
339    let q = rnpb[1][0] * xs + rnpb[1][1] * ys + rnpb[1][2] * zs;
340    let eo = if p != 0.0 || q != 0.0 {
341        s - q.atan2(p)
342    } else {
343        s
344    };
345    Ok(eo)
346}
erfa_rust/G11_safe.rs

erfa_rust/
G11_safe.rs
