erfa_rust/
G29_safe.rs

1// G29_safe  CIO locator & spherical/Cartesian
2//   s00.c   → eraS00_safe
3//   s00a.c  → eraS00a_safe
4//   s00b.c  → eraS00b_safe
5//   s2c.c   → eraS2c_safe
6//   s2p.c   → eraS2p_safe
7//   s2pv.c  → eraS2pv_safe
8//   s2xpv.c → eraS2xpv_safe
9
10use crate::G15_safe::{
11    eraFad03_safe, eraFae03_safe, eraFaf03_safe, eraFal03_safe, eraFalp03_safe, eraFaom03_safe,
12    eraFapa03_safe, eraFave03_safe,
13};
14use crate::G26_safe::{eraPnm00a_safe, eraPnm00b_safe};
15use crate::G6_safe::eraBpn2xy_safe;
16use crate::H1_safe::{ERFA_DAS2R, ERFA_DJ00, ERFA_DJC};
17
18pub type ErfaResult<T> = Result<T, ()>;
19
20// ===================================
21// G29/s00.c
22// ===================================
23// Compute CIO locator s (IAU 2000), given date and X,Y.
24pub fn eraS00_safe(date1: f64, date2: f64, x: f64, y: f64) -> ErfaResult<f64> {
25    // Time since J2000.0 (Julian centuries)
26    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
27
28    // Series for s + X*Y/2
29    #[derive(Clone, Copy)]
30    struct Term {
31        nfa: [i32; 8], // coeffs of l,l',F,D,Om,LVe,LE,pA
32        s: f64,        // sine coefficient (µas)
33        c: f64,        // cosine coefficient (µas)
34    }
35
36    // Polynomial coefficients (µas)
37    const SP: [f64; 6] = [
38        94.00e-6,
39        3808.35e-6,
40        -119.94e-6,
41        -72574.09e-6,
42        27.70e-6,
43        15.61e-6,
44    ];
45
46    // Terms order t^0
47    const S0: [Term; 33] = [
48        Term {
49            nfa: [0, 0, 0, 0, 1, 0, 0, 0],
50            s: -2640.73e-6,
51            c: 0.39e-6,
52        },
53        Term {
54            nfa: [0, 0, 0, 0, 2, 0, 0, 0],
55            s: -63.53e-6,
56            c: 0.02e-6,
57        },
58        Term {
59            nfa: [0, 0, 2, -2, 3, 0, 0, 0],
60            s: -11.75e-6,
61            c: -0.01e-6,
62        },
63        Term {
64            nfa: [0, 0, 2, -2, 1, 0, 0, 0],
65            s: -11.21e-6,
66            c: -0.01e-6,
67        },
68        Term {
69            nfa: [0, 0, 2, -2, 2, 0, 0, 0],
70            s: 4.57e-6,
71            c: 0.00e-6,
72        },
73        Term {
74            nfa: [0, 0, 2, 0, 3, 0, 0, 0],
75            s: -2.02e-6,
76            c: 0.00e-6,
77        },
78        Term {
79            nfa: [0, 0, 2, 0, 1, 0, 0, 0],
80            s: -1.98e-6,
81            c: 0.00e-6,
82        },
83        Term {
84            nfa: [0, 0, 0, 0, 3, 0, 0, 0],
85            s: 1.72e-6,
86            c: 0.00e-6,
87        },
88        Term {
89            nfa: [0, 1, 0, 0, 1, 0, 0, 0],
90            s: 1.41e-6,
91            c: 0.01e-6,
92        },
93        Term {
94            nfa: [0, 1, 0, 0, -1, 0, 0, 0],
95            s: 1.26e-6,
96            c: 0.01e-6,
97        },
98        Term {
99            nfa: [1, 0, 0, 0, -1, 0, 0, 0],
100            s: 0.63e-6,
101            c: 0.00e-6,
102        },
103        Term {
104            nfa: [1, 0, 0, 0, 1, 0, 0, 0],
105            s: 0.63e-6,
106            c: 0.00e-6,
107        },
108        Term {
109            nfa: [0, 1, 2, -2, 3, 0, 0, 0],
110            s: -0.46e-6,
111            c: 0.00e-6,
112        },
113        Term {
114            nfa: [0, 1, 2, -2, 1, 0, 0, 0],
115            s: -0.45e-6,
116            c: 0.00e-6,
117        },
118        Term {
119            nfa: [0, 0, 4, -4, 4, 0, 0, 0],
120            s: -0.36e-6,
121            c: 0.00e-6,
122        },
123        Term {
124            nfa: [0, 0, 1, -1, 1, -8, 12, 0],
125            s: 0.24e-6,
126            c: 0.12e-6,
127        },
128        Term {
129            nfa: [0, 0, 2, 0, 0, 0, 0, 0],
130            s: -0.32e-6,
131            c: 0.00e-6,
132        },
133        Term {
134            nfa: [0, 0, 2, 0, 2, 0, 0, 0],
135            s: -0.28e-6,
136            c: 0.00e-6,
137        },
138        Term {
139            nfa: [1, 0, 2, 0, 3, 0, 0, 0],
140            s: -0.27e-6,
141            c: 0.00e-6,
142        },
143        Term {
144            nfa: [1, 0, 2, 0, 1, 0, 0, 0],
145            s: -0.26e-6,
146            c: 0.00e-6,
147        },
148        Term {
149            nfa: [0, 0, 2, -2, 0, 0, 0, 0],
150            s: 0.21e-6,
151            c: 0.00e-6,
152        },
153        Term {
154            nfa: [0, 1, -2, 2, -3, 0, 0, 0],
155            s: -0.19e-6,
156            c: 0.00e-6,
157        },
158        Term {
159            nfa: [0, 1, -2, 2, -1, 0, 0, 0],
160            s: -0.18e-6,
161            c: 0.00e-6,
162        },
163        Term {
164            nfa: [0, 0, 0, 0, 0, 8, -13, -1],
165            s: 0.10e-6,
166            c: -0.05e-6,
167        },
168        Term {
169            nfa: [0, 0, 0, 2, 0, 0, 0, 0],
170            s: -0.15e-6,
171            c: 0.00e-6,
172        },
173        Term {
174            nfa: [2, 0, -2, 0, -1, 0, 0, 0],
175            s: 0.14e-6,
176            c: 0.00e-6,
177        },
178        Term {
179            nfa: [0, 1, 2, -2, 2, 0, 0, 0],
180            s: 0.14e-6,
181            c: 0.00e-6,
182        },
183        Term {
184            nfa: [1, 0, 0, -2, 1, 0, 0, 0],
185            s: -0.14e-6,
186            c: 0.00e-6,
187        },
188        Term {
189            nfa: [1, 0, 0, -2, -1, 0, 0, 0],
190            s: -0.14e-6,
191            c: 0.00e-6,
192        },
193        Term {
194            nfa: [0, 0, 4, -2, 4, 0, 0, 0],
195            s: -0.13e-6,
196            c: 0.00e-6,
197        },
198        Term {
199            nfa: [0, 0, 2, -2, 4, 0, 0, 0],
200            s: 0.11e-6,
201            c: 0.00e-6,
202        },
203        Term {
204            nfa: [1, 0, -2, 0, -3, 0, 0, 0],
205            s: -0.11e-6,
206            c: 0.00e-6,
207        },
208        Term {
209            nfa: [1, 0, -2, 0, -1, 0, 0, 0],
210            s: -0.11e-6,
211            c: 0.00e-6,
212        },
213    ];
214
215    // Terms order t^1
216    const S1: [Term; 3] = [
217        Term {
218            nfa: [0, 0, 0, 0, 2, 0, 0, 0],
219            s: -0.07e-6,
220            c: 3.57e-6,
221        },
222        Term {
223            nfa: [0, 0, 0, 0, 1, 0, 0, 0],
224            s: 1.71e-6,
225            c: -0.03e-6,
226        },
227        Term {
228            nfa: [0, 0, 2, -2, 3, 0, 0, 0],
229            s: 0.00e-6,
230            c: 0.48e-6,
231        },
232    ];
233
234    // Terms order t^2
235    const S2: [Term; 25] = [
236        Term {
237            nfa: [0, 0, 0, 0, 1, 0, 0, 0],
238            s: 743.53e-6,
239            c: -0.17e-6,
240        },
241        Term {
242            nfa: [0, 0, 2, -2, 2, 0, 0, 0],
243            s: 56.91e-6,
244            c: 0.06e-6,
245        },
246        Term {
247            nfa: [0, 0, 2, 0, 2, 0, 0, 0],
248            s: 9.84e-6,
249            c: -0.01e-6,
250        },
251        Term {
252            nfa: [0, 0, 0, 0, 2, 0, 0, 0],
253            s: -8.85e-6,
254            c: 0.01e-6,
255        },
256        Term {
257            nfa: [0, 1, 0, 0, 0, 0, 0, 0],
258            s: -6.38e-6,
259            c: -0.05e-6,
260        },
261        Term {
262            nfa: [1, 0, 0, 0, 0, 0, 0, 0],
263            s: -3.07e-6,
264            c: 0.00e-6,
265        },
266        Term {
267            nfa: [0, 1, 2, -2, 2, 0, 0, 0],
268            s: 2.23e-6,
269            c: 0.00e-6,
270        },
271        Term {
272            nfa: [0, 0, 2, 0, 1, 0, 0, 0],
273            s: 1.67e-6,
274            c: 0.00e-6,
275        },
276        Term {
277            nfa: [1, 0, 2, 0, 2, 0, 0, 0],
278            s: 1.30e-6,
279            c: 0.00e-6,
280        },
281        Term {
282            nfa: [0, 1, -2, 2, -2, 0, 0, 0],
283            s: 0.93e-6,
284            c: 0.00e-6,
285        },
286        Term {
287            nfa: [1, 0, 0, -2, 0, 0, 0, 0],
288            s: 0.68e-6,
289            c: 0.00e-6,
290        },
291        Term {
292            nfa: [0, 0, 2, -2, 1, 0, 0, 0],
293            s: -0.55e-6,
294            c: 0.00e-6,
295        },
296        Term {
297            nfa: [1, 0, -2, 0, -2, 0, 0, 0],
298            s: 0.53e-6,
299            c: 0.00e-6,
300        },
301        Term {
302            nfa: [0, 0, 0, 2, 0, 0, 0, 0],
303            s: -0.27e-6,
304            c: 0.00e-6,
305        },
306        Term {
307            nfa: [1, 0, 0, 0, 1, 0, 0, 0],
308            s: -0.27e-6,
309            c: 0.00e-6,
310        },
311        Term {
312            nfa: [1, 0, -2, -2, -2, 0, 0, 0],
313            s: -0.26e-6,
314            c: 0.00e-6,
315        },
316        Term {
317            nfa: [1, 0, 0, 0, -1, 0, 0, 0],
318            s: -0.25e-6,
319            c: 0.00e-6,
320        },
321        Term {
322            nfa: [1, 0, 2, 0, 1, 0, 0, 0],
323            s: 0.22e-6,
324            c: 0.00e-6,
325        },
326        Term {
327            nfa: [2, 0, 0, -2, 0, 0, 0, 0],
328            s: -0.21e-6,
329            c: 0.00e-6,
330        },
331        Term {
332            nfa: [2, 0, -2, 0, -1, 0, 0, 0],
333            s: 0.20e-6,
334            c: 0.00e-6,
335        },
336        Term {
337            nfa: [0, 0, 2, 2, 2, 0, 0, 0],
338            s: 0.17e-6,
339            c: 0.00e-6,
340        },
341        Term {
342            nfa: [2, 0, 2, 0, 2, 0, 0, 0],
343            s: 0.13e-6,
344            c: 0.00e-6,
345        },
346        Term {
347            nfa: [2, 0, 0, 0, 0, 0, 0, 0],
348            s: -0.13e-6,
349            c: 0.00e-6,
350        },
351        Term {
352            nfa: [1, 0, 2, -2, 2, 0, 0, 0],
353            s: -0.12e-6,
354            c: 0.00e-6,
355        },
356        Term {
357            nfa: [0, 0, 2, 0, 0, 0, 0, 0],
358            s: -0.11e-6,
359            c: 0.00e-6,
360        },
361    ];
362
363    // Terms order t^3
364    const S3: [Term; 4] = [
365        Term {
366            nfa: [0, 0, 0, 0, 1, 0, 0, 0],
367            s: 0.30e-6,
368            c: -23.51e-6,
369        },
370        Term {
371            nfa: [0, 0, 2, -2, 2, 0, 0, 0],
372            s: -0.03e-6,
373            c: -1.39e-6,
374        },
375        Term {
376            nfa: [0, 0, 2, 0, 2, 0, 0, 0],
377            s: -0.01e-6,
378            c: -0.24e-6,
379        },
380        Term {
381            nfa: [0, 0, 0, 0, 2, 0, 0, 0],
382            s: 0.00e-6,
383            c: 0.22e-6,
384        },
385    ];
386
387    // Terms order t^4
388    const S4: [Term; 1] = [Term {
389        nfa: [0, 0, 0, 0, 1, 0, 0, 0],
390        s: -0.26e-6,
391        c: -0.01e-6,
392    }];
393
394    // Fundamental arguments (radians)
395    let mut fa = [0.0_f64; 8];
396    fa[0] = eraFal03_safe(t)?;
397    fa[1] = eraFalp03_safe(t)?;
398    fa[2] = eraFaf03_safe(t)?;
399    fa[3] = eraFad03_safe(t)?;
400    fa[4] = eraFaom03_safe(t)?;
401    fa[5] = eraFave03_safe(t)?;
402    fa[6] = eraFae03_safe(t)?;
403    fa[7] = eraFapa03_safe(t)?;
404
405    // Evaluate series
406    let mut w = [SP[0], SP[1], SP[2], SP[3], SP[4], SP[5]];
407
408    #[inline(always)]
409    fn accum(arr: &[Term], idx: usize, fa: &[f64; 8], w: &mut [f64; 6]) {
410        for term in arr {
411            let mut a = 0.0_f64;
412            for j in 0..8 {
413                a += (term.nfa[j] as f64) * fa[j];
414            }
415            w[idx] += term.s * a.sin() + term.c * a.cos();
416        }
417    }
418    accum(&S0, 0, &fa, &mut w);
419    accum(&S1, 1, &fa, &mut w);
420    accum(&S2, 2, &fa, &mut w);
421    accum(&S3, 3, &fa, &mut w);
422    accum(&S4, 4, &fa, &mut w);
423
424    // Final value (radians)
425    let s = (w[0] + (w[1] + (w[2] + (w[3] + (w[4] + w[5] * t) * t) * t) * t) * t) * ERFA_DAS2R
426        - x * y / 2.0;
427
428    Ok(s)
429}
430
431// ===================================
432// G29/s00a.c
433// ===================================
434// CIO locator s using full IAU 2000A model.
435pub fn eraS00a_safe(date1: f64, date2: f64) -> ErfaResult<f64> {
436    let rbpn = eraPnm00a_safe(date1, date2)?;
437    let (x, y) = eraBpn2xy_safe(&rbpn)?;
438    eraS00_safe(date1, date2, x, y)
439}
440
441// ===================================
442// G29/s00b.c
443// ===================================
444// CIO locator s using truncated IAU 2000B model.
445pub fn eraS00b_safe(date1: f64, date2: f64) -> ErfaResult<f64> {
446    let rbpn = eraPnm00b_safe(date1, date2)?;
447    let (x, y) = eraBpn2xy_safe(&rbpn)?;
448    eraS00_safe(date1, date2, x, y)
449}
450
451// ===================================
452// G29/s2c.c
453// ===================================
454// Convert spherical angles to Cartesian unit vector.
455pub fn eraS2c_safe(theta: f64, phi: f64) -> ErfaResult<[f64; 3]> {
456    let cp = phi.cos();
457    Ok([theta.cos() * cp, theta.sin() * cp, phi.sin()])
458}
459
460// ===================================
461// G29/s2p.c
462// ===================================
463// Convert spherical polar to position vector with radius r.
464pub fn eraS2p_safe(theta: f64, phi: f64, r: f64) -> ErfaResult<[f64; 3]> {
465    let u = eraS2c_safe(theta, phi)?;
466    Ok([r * u[0], r * u[1], r * u[2]])
467}
468
469// ===================================
470// G29/s2pv.c
471// ===================================
472// Convert spherical angles and rates to pv-vector.
473pub fn eraS2pv_safe(
474    theta: f64,
475    phi: f64,
476    r: f64,
477    td: f64,
478    pd: f64,
479    rd: f64,
480) -> ErfaResult<[[f64; 3]; 2]> {
481    let st = theta.sin();
482    let ct = theta.cos();
483    let sp = phi.sin();
484    let cp = phi.cos();
485
486    let rcp = r * cp;
487    let x = rcp * ct;
488    let y = rcp * st;
489    let rpd = r * pd;
490    let w = rpd * sp - cp * rd;
491
492    let p = [x, y, r * sp];
493    let v = [-y * td - w * ct, x * td - w * st, rpd * cp + sp * rd];
494
495    Ok([p, v])
496}
497
498// ===================================
499// G29/s2xpv.c
500// ===================================
501// Scale pv-vector components by two scalars (position by s1, velocity by s2).
502pub fn eraS2xpv_safe(s1: f64, s2: f64, pv: &[[f64; 3]; 2]) -> ErfaResult<[[f64; 3]; 2]> {
503    let p = [s1 * pv[0][0], s1 * pv[0][1], s1 * pv[0][2]];
504    let v = [s2 * pv[1][0], s2 * pv[1][1], s2 * pv[1][2]];
505    Ok([p, v])
506}
erfa_rust/G29_safe.rs

erfa_rust/
G29_safe.rs
