erfa_rust/
G30_safe.rs

1// G30
2//   s06.c    → eraS06_safe
3//   s06a.c   → eraS06a_safe
4//   sepp.c   → eraSepp_safe
5//   seps.c   → eraSeps_safe
6//   sp00.c   → eraSp00_safe
7//   starpm.c → eraStarpm_safe
8//   starpv.c → eraStarpv_safe
9//   sxp.c    → eraSxp_safe
10//   sxpv.c   → eraSxpv_safe
11
12use crate::G15_safe::{
13    eraFad03_safe, eraFae03_safe, eraFaf03_safe, eraFal03_safe, eraFalp03_safe, eraFaom03_safe,
14    eraFapa03_safe, eraFave03_safe,
15};
16use crate::G24_safe::{eraPdp_safe, eraPm_safe};
17use crate::G25_safe::{eraPmp_safe, eraPn_safe};
18use crate::G26_safe::{eraPnm06a_safe, eraPpp_safe};
19use crate::G27_safe::{eraPvstar_safe, eraPvu_safe, eraPxp_safe};
20use crate::G29_safe::{eraS2c_safe, eraS2pv_safe, eraS2xpv_safe};
21use crate::G35_safe::eraZp_safe;
22use crate::G6_safe::eraBpn2xy_safe;
23use crate::H1_safe::{
24    ERFA_DAS2R, ERFA_DAU, ERFA_DAYSEC, ERFA_DC, ERFA_DJ00, ERFA_DJC, ERFA_DJY, ERFA_DR2AS,
25};
26
27pub type ErfaResult<T> = Result<T, ()>;
28
29//----------------------------------------------------------------------
30// G30/s06.c
31//----------------------------------------------------------------------
32// Internal data for s06 series
33#[derive(Clone, Copy)]
34struct TERM {
35    nfa: [i32; 8], // coefficients of l,l',F,D,Om,LVe,LE,pA
36    s: f64,        // sine coefficient (µas)
37    c: f64,        // cosine coefficient (µas)
38}
39
40const SP: [f64; 6] = [
41    94.00e-6,
42    3808.65e-6,
43    -122.68e-6,
44    -72574.11e-6,
45    27.98e-6,
46    15.62e-6,
47];
48
49const S0: [TERM; 33] = [
50    TERM {
51        nfa: [0, 0, 0, 0, 1, 0, 0, 0],
52        s: -2640.73e-6,
53        c: 0.39e-6,
54    },
55    TERM {
56        nfa: [0, 0, 0, 0, 2, 0, 0, 0],
57        s: -63.53e-6,
58        c: 0.02e-6,
59    },
60    TERM {
61        nfa: [0, 0, 2, -2, 3, 0, 0, 0],
62        s: -11.75e-6,
63        c: -0.01e-6,
64    },
65    TERM {
66        nfa: [0, 0, 2, -2, 1, 0, 0, 0],
67        s: -11.21e-6,
68        c: -0.01e-6,
69    },
70    TERM {
71        nfa: [0, 0, 2, -2, 2, 0, 0, 0],
72        s: 4.57e-6,
73        c: 0.00e-6,
74    },
75    TERM {
76        nfa: [0, 0, 2, 0, 3, 0, 0, 0],
77        s: -2.02e-6,
78        c: 0.00e-6,
79    },
80    TERM {
81        nfa: [0, 0, 2, 0, 1, 0, 0, 0],
82        s: -1.98e-6,
83        c: 0.00e-6,
84    },
85    TERM {
86        nfa: [0, 0, 0, 0, 3, 0, 0, 0],
87        s: 1.72e-6,
88        c: 0.00e-6,
89    },
90    TERM {
91        nfa: [0, 1, 0, 0, 1, 0, 0, 0],
92        s: 1.41e-6,
93        c: 0.01e-6,
94    },
95    TERM {
96        nfa: [0, 1, 0, 0, -1, 0, 0, 0],
97        s: 1.26e-6,
98        c: 0.01e-6,
99    },
100    // 11-20
101    TERM {
102        nfa: [1, 0, 0, 0, -1, 0, 0, 0],
103        s: 0.63e-6,
104        c: 0.00e-6,
105    },
106    TERM {
107        nfa: [1, 0, 0, 0, 1, 0, 0, 0],
108        s: 0.63e-6,
109        c: 0.00e-6,
110    },
111    TERM {
112        nfa: [0, 1, 2, -2, 3, 0, 0, 0],
113        s: -0.46e-6,
114        c: 0.00e-6,
115    },
116    TERM {
117        nfa: [0, 1, 2, -2, 1, 0, 0, 0],
118        s: -0.45e-6,
119        c: 0.00e-6,
120    },
121    TERM {
122        nfa: [0, 0, 4, -4, 4, 0, 0, 0],
123        s: -0.36e-6,
124        c: 0.00e-6,
125    },
126    TERM {
127        nfa: [0, 0, 1, -1, 1, -8, 12, 0],
128        s: 0.24e-6,
129        c: 0.12e-6,
130    },
131    TERM {
132        nfa: [0, 0, 2, 0, 0, 0, 0, 0],
133        s: -0.32e-6,
134        c: 0.00e-6,
135    },
136    TERM {
137        nfa: [0, 0, 2, 0, 2, 0, 0, 0],
138        s: -0.28e-6,
139        c: 0.00e-6,
140    },
141    TERM {
142        nfa: [1, 0, 2, 0, 3, 0, 0, 0],
143        s: -0.27e-6,
144        c: 0.00e-6,
145    },
146    TERM {
147        nfa: [1, 0, 2, 0, 1, 0, 0, 0],
148        s: -0.26e-6,
149        c: 0.00e-6,
150    },
151    // 21-30
152    TERM {
153        nfa: [0, 0, 2, -2, 0, 0, 0, 0],
154        s: 0.21e-6,
155        c: 0.00e-6,
156    },
157    TERM {
158        nfa: [0, 1, -2, 2, -3, 0, 0, 0],
159        s: -0.19e-6,
160        c: 0.00e-6,
161    },
162    TERM {
163        nfa: [0, 1, -2, 2, -1, 0, 0, 0],
164        s: -0.18e-6,
165        c: 0.00e-6,
166    },
167    TERM {
168        nfa: [0, 0, 0, 0, 0, 8, -13, -1],
169        s: 0.10e-6,
170        c: -0.05e-6,
171    },
172    TERM {
173        nfa: [0, 0, 0, 2, 0, 0, 0, 0],
174        s: -0.15e-6,
175        c: 0.00e-6,
176    },
177    TERM {
178        nfa: [2, 0, -2, 0, -1, 0, 0, 0],
179        s: 0.14e-6,
180        c: 0.00e-6,
181    },
182    TERM {
183        nfa: [0, 1, 2, -2, 2, 0, 0, 0],
184        s: 0.14e-6,
185        c: 0.00e-6,
186    },
187    TERM {
188        nfa: [1, 0, 0, -2, 1, 0, 0, 0],
189        s: -0.14e-6,
190        c: 0.00e-6,
191    },
192    TERM {
193        nfa: [1, 0, 0, -2, -1, 0, 0, 0],
194        s: -0.14e-6,
195        c: 0.00e-6,
196    },
197    TERM {
198        nfa: [0, 0, 4, -2, 4, 0, 0, 0],
199        s: -0.13e-6,
200        c: 0.00e-6,
201    },
202    // 31-33
203    TERM {
204        nfa: [0, 0, 2, -2, 4, 0, 0, 0],
205        s: 0.11e-6,
206        c: 0.00e-6,
207    },
208    TERM {
209        nfa: [1, 0, -2, 0, -3, 0, 0, 0],
210        s: -0.11e-6,
211        c: 0.00e-6,
212    },
213    TERM {
214        nfa: [1, 0, -2, 0, -1, 0, 0, 0],
215        s: -0.11e-6,
216        c: 0.00e-6,
217    },
218];
219
220// t^1
221const S1: [TERM; 3] = [
222    TERM {
223        nfa: [0, 0, 0, 0, 2, 0, 0, 0],
224        s: -0.07e-6,
225        c: 3.57e-6,
226    },
227    TERM {
228        nfa: [0, 0, 0, 0, 1, 0, 0, 0],
229        s: 1.73e-6,
230        c: -0.03e-6,
231    },
232    TERM {
233        nfa: [0, 0, 2, -2, 3, 0, 0, 0],
234        s: 0.00e-6,
235        c: 0.48e-6,
236    },
237];
238
239// t^2
240const S2: [TERM; 25] = [
241    // 1-10
242    TERM {
243        nfa: [0, 0, 0, 0, 1, 0, 0, 0],
244        s: 743.52e-6,
245        c: -0.17e-6,
246    },
247    TERM {
248        nfa: [0, 0, 2, -2, 2, 0, 0, 0],
249        s: 56.91e-6,
250        c: 0.06e-6,
251    },
252    TERM {
253        nfa: [0, 0, 2, 0, 2, 0, 0, 0],
254        s: 9.84e-6,
255        c: -0.01e-6,
256    },
257    TERM {
258        nfa: [0, 0, 0, 0, 2, 0, 0, 0],
259        s: -8.85e-6,
260        c: 0.01e-6,
261    },
262    TERM {
263        nfa: [0, 1, 0, 0, 0, 0, 0, 0],
264        s: -6.38e-6,
265        c: -0.05e-6,
266    },
267    TERM {
268        nfa: [1, 0, 0, 0, 0, 0, 0, 0],
269        s: -3.07e-6,
270        c: 0.00e-6,
271    },
272    TERM {
273        nfa: [0, 1, 2, -2, 2, 0, 0, 0],
274        s: 2.23e-6,
275        c: 0.00e-6,
276    },
277    TERM {
278        nfa: [0, 0, 2, 0, 1, 0, 0, 0],
279        s: 1.67e-6,
280        c: 0.00e-6,
281    },
282    TERM {
283        nfa: [1, 0, 2, 0, 2, 0, 0, 0],
284        s: 1.30e-6,
285        c: 0.00e-6,
286    },
287    TERM {
288        nfa: [0, 1, -2, 2, -2, 0, 0, 0],
289        s: 0.93e-6,
290        c: 0.00e-6,
291    },
292    // 11-20
293    TERM {
294        nfa: [1, 0, 0, -2, 0, 0, 0, 0],
295        s: 0.68e-6,
296        c: 0.00e-6,
297    },
298    TERM {
299        nfa: [0, 0, 2, -2, 1, 0, 0, 0],
300        s: -0.55e-6,
301        c: 0.00e-6,
302    },
303    TERM {
304        nfa: [1, 0, -2, 0, -2, 0, 0, 0],
305        s: 0.53e-6,
306        c: 0.00e-6,
307    },
308    TERM {
309        nfa: [0, 0, 0, 2, 0, 0, 0, 0],
310        s: -0.27e-6,
311        c: 0.00e-6,
312    },
313    TERM {
314        nfa: [1, 0, 0, 0, 1, 0, 0, 0],
315        s: -0.27e-6,
316        c: 0.00e-6,
317    },
318    TERM {
319        nfa: [1, 0, -2, -2, -2, 0, 0, 0],
320        s: -0.26e-6,
321        c: 0.00e-6,
322    },
323    TERM {
324        nfa: [1, 0, 0, 0, -1, 0, 0, 0],
325        s: -0.25e-6,
326        c: 0.00e-6,
327    },
328    TERM {
329        nfa: [1, 0, 2, 0, 1, 0, 0, 0],
330        s: 0.22e-6,
331        c: 0.00e-6,
332    },
333    TERM {
334        nfa: [2, 0, 0, -2, 0, 0, 0, 0],
335        s: -0.21e-6,
336        c: 0.00e-6,
337    },
338    TERM {
339        nfa: [2, 0, -2, 0, -1, 0, 0, 0],
340        s: 0.20e-6,
341        c: 0.00e-6,
342    },
343    // 21-25
344    TERM {
345        nfa: [0, 0, 2, 2, 2, 0, 0, 0],
346        s: 0.17e-6,
347        c: 0.00e-6,
348    },
349    TERM {
350        nfa: [2, 0, 2, 0, 2, 0, 0, 0],
351        s: 0.13e-6,
352        c: 0.00e-6,
353    },
354    TERM {
355        nfa: [2, 0, 0, 0, 0, 0, 0, 0],
356        s: -0.13e-6,
357        c: 0.00e-6,
358    },
359    TERM {
360        nfa: [1, 0, 2, -2, 2, 0, 0, 0],
361        s: -0.12e-6,
362        c: 0.00e-6,
363    },
364    TERM {
365        nfa: [0, 0, 2, 0, 0, 0, 0, 0],
366        s: -0.11e-6,
367        c: 0.00e-6,
368    },
369];
370
371// t^3
372const S3: [TERM; 4] = [
373    TERM {
374        nfa: [0, 0, 0, 0, 1, 0, 0, 0],
375        s: 0.30e-6,
376        c: -23.42e-6,
377    },
378    TERM {
379        nfa: [0, 0, 2, -2, 2, 0, 0, 0],
380        s: -0.03e-6,
381        c: -1.46e-6,
382    },
383    TERM {
384        nfa: [0, 0, 2, 0, 2, 0, 0, 0],
385        s: -0.01e-6,
386        c: -0.25e-6,
387    },
388    TERM {
389        nfa: [0, 0, 0, 0, 2, 0, 0, 0],
390        s: 0.00e-6,
391        c: 0.23e-6,
392    },
393];
394
395// t^4
396const S4: [TERM; 1] = [TERM {
397    nfa: [0, 0, 0, 0, 1, 0, 0, 0],
398    s: -0.26e-6,
399    c: -0.01e-6,
400}];
401
402// CIO locator s (IAU 2006/2000A), given CIP X,Y.
403pub fn eraS06_safe(date1: f64, date2: f64, x: f64, y: f64) -> ErfaResult<f64> {
404    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
405
406    let mut fa = [0.0_f64; 8];
407    fa[0] = eraFal03_safe(t)?; // l
408    fa[1] = eraFalp03_safe(t)?; // l'
409    fa[2] = eraFaf03_safe(t)?; // F
410    fa[3] = eraFad03_safe(t)?; // D
411    fa[4] = eraFaom03_safe(t)?; // Om
412    fa[5] = eraFave03_safe(t)?; // LVe
413    fa[6] = eraFae03_safe(t)?; // LE
414    fa[7] = eraFapa03_safe(t)?; // pA
415
416    // evaluate series for s + X*Y/2
417    let mut w0 = SP[0];
418    let mut w1 = SP[1];
419    let mut w2 = SP[2];
420    let mut w3 = SP[3];
421    let mut w4 = SP[4];
422    let w5 = SP[5];
423
424    for term in S0.iter().rev() {
425        let mut a = 0.0_f64;
426        for j in 0..8 {
427            a += term.nfa[j] as f64 * fa[j];
428        }
429        w0 += term.s * a.sin() + term.c * a.cos();
430    }
431
432    for term in S1.iter().rev() {
433        let mut a = 0.0_f64;
434        for j in 0..8 {
435            a += term.nfa[j] as f64 * fa[j];
436        }
437        w1 += term.s * a.sin() + term.c * a.cos();
438    }
439
440    for term in S2.iter().rev() {
441        let mut a = 0.0_f64;
442        for j in 0..8 {
443            a += term.nfa[j] as f64 * fa[j];
444        }
445        w2 += term.s * a.sin() + term.c * a.cos();
446    }
447
448    for term in S3.iter().rev() {
449        let mut a = 0.0_f64;
450        for j in 0..8 {
451            a += term.nfa[j] as f64 * fa[j];
452        }
453        w3 += term.s * a.sin() + term.c * a.cos();
454    }
455
456    for term in S4.iter().rev() {
457        let mut a = 0.0_f64;
458        for j in 0..8 {
459            a += term.nfa[j] as f64 * fa[j];
460        }
461        w4 += term.s * a.sin() + term.c * a.cos();
462    }
463
464    let s = (w0 + (w1 + (w2 + (w3 + (w4 + w5 * t) * t) * t) * t) * t) * ERFA_DAS2R - x * y / 2.0;
465
466    Ok(s)
467}
468
469//----------------------------------------------------------------------
470// G30/s06a.c
471//----------------------------------------------------------------------
472// CIO locator s, IAU 2006/2000A (auto X,Y).
473pub fn eraS06a_safe(date1: f64, date2: f64) -> ErfaResult<f64> {
474    let rnpb = eraPnm06a_safe(date1, date2)?;
475    let (x, y) = eraBpn2xy_safe(&rnpb)?;
476    eraS06_safe(date1, date2, x, y)
477}
478
479//----------------------------------------------------------------------
480// G30/sepp.c
481//----------------------------------------------------------------------
482// Angular separation between two p-vectors.
483pub fn eraSepp_safe(a: &[f64; 3], b: &[f64; 3]) -> ErfaResult<f64> {
484    let axb = eraPxp_safe(a, b)?;
485    let ss = eraPm_safe(&axb)?;
486    let cs = eraPdp_safe(a, b)?;
487    let ang = if ss != 0.0 || cs != 0.0 {
488        ss.atan2(cs)
489    } else {
490        0.0
491    };
492    Ok(ang)
493}
494
495//----------------------------------------------------------------------
496// G30/seps.c
497//----------------------------------------------------------------------
498// Angular separation between two spherical positions.
499pub fn eraSeps_safe(al: f64, ap: f64, bl: f64, bp: f64) -> ErfaResult<f64> {
500    let ac = eraS2c_safe(al, ap)?;
501    let bc = eraS2c_safe(bl, bp)?;
502    eraSepp_safe(&ac, &bc)
503}
504
505//----------------------------------------------------------------------
506// G30/sp00.c
507//----------------------------------------------------------------------
508// TIO locator s′ (linear model).
509pub fn eraSp00_safe(date1: f64, date2: f64) -> ErfaResult<f64> {
510    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
511    Ok(-47e-6 * t * ERFA_DAS2R)
512}
513
514//----------------------------------------------------------------------
515// G30/starpm.c
516//----------------------------------------------------------------------
517// Propagate star catalog data for space motion.
518pub fn eraStarpm_safe(
519    ra1: f64,
520    dec1: f64,
521    pmr1: f64,
522    pmd1: f64,
523    px1: f64,
524    rv1: f64,
525    ep1a: f64,
526    ep1b: f64,
527    ep2a: f64,
528    ep2b: f64,
529) -> ErfaResult<((f64, f64, f64, f64, f64, f64), i32)> {
530    let (pv1, j1) = eraStarpv_safe(ra1, dec1, pmr1, pmd1, px1, rv1)?;
531    let tl1 = eraPm_safe(&pv1[0])? / ERFA_DC;
532    let dt = (ep2a - ep1a) + (ep2b - ep1b);
533    let pv = eraPvu_safe(dt + tl1, &pv1)?;
534
535    let r2 = eraPdp_safe(&pv[0], &pv[0])?;
536    let rdv = eraPdp_safe(&pv[0], &pv[1])?;
537    let v2 = eraPdp_safe(&pv[1], &pv[1])?;
538    let c2mv2 = ERFA_DC * ERFA_DC - v2;
539    if c2mv2 <= 0.0 {
540        return Ok(((0.0, 0.0, 0.0, 0.0, 0.0, 0.0), -1));
541    }
542    let tl2 = (-rdv + (rdv * rdv + c2mv2 * r2).sqrt()) / c2mv2;
543
544    let pv2 = eraPvu_safe(dt + (tl1 - tl2), &pv1)?;
545    let ((ra2, dec2, pmr2, pmd2, px2, rv2), j2) = eraPvstar_safe(&pv2)?;
546    let j = if j2 == 0 { j1 } else { -1 };
547    Ok(((ra2, dec2, pmr2, pmd2, px2, rv2), j))
548}
549
550//----------------------------------------------------------------------
551// G30/starpv.c
552//----------------------------------------------------------------------
553// Catalog data to position/velocity vector.
554pub fn eraStarpv_safe(
555    ra: f64,
556    dec: f64,
557    pmr: f64,
558    pmd: f64,
559    px: f64,
560    rv: f64,
561) -> ErfaResult<([[f64; 3]; 2], i32)> {
562    const PXMIN: f64 = 1e-7;
563    const VMAX: f64 = 0.5;
564    const IMAX: i32 = 100;
565
566    let mut warn = 0;
567
568    let w = if px >= PXMIN {
569        px
570    } else {
571        warn += 1;
572        PXMIN
573    };
574    let r = ERFA_DR2AS / w;
575
576    let rd = ERFA_DAYSEC * rv * 1e3 / ERFA_DAU;
577    let rad = pmr / ERFA_DJY;
578    let decd = pmd / ERFA_DJY;
579
580    let pv_init = eraS2pv_safe(ra, dec, r, rad, decd, rd)?;
581
582    let mut pv = pv_init;
583    let v = eraPm_safe(&pv[1])?;
584    if v / ERFA_DC > VMAX {
585        pv[1] = eraZp_safe();
586        warn += 2;
587    }
588
589    let (_r_u, pu) = eraPn_safe(&pv[0])?;
590    let vsr = eraPdp_safe(&pu, &pv[1])?;
591    let usr = eraSxp_safe(vsr, &pu)?;
592
593    let ust = eraPmp_safe(&pv[1], &usr)?;
594    let vst = eraPm_safe(&ust)?;
595
596    let betsr = vsr / ERFA_DC;
597    let betst = vst / ERFA_DC;
598
599    let mut betr = betsr;
600    let mut bett = betst;
601    let mut d = 1.0_f64;
602    let mut del = 0.0_f64;
603    let mut od = 0.0_f64;
604    let mut odel = 0.0_f64;
605    let mut odd = 0.0_f64;
606    let mut oddel = 0.0_f64;
607
608    let mut i = 0;
609    while i < IMAX {
610        d = 1.0 + betr;
611        let w2 = betr * betr + bett * bett;
612        del = -w2 / ((1.0 - w2).sqrt() + 1.0);
613        betr = d * betsr + del;
614        bett = d * betst;
615
616        if i > 0 {
617            let dd = (d - od).abs();
618            let ddel = (del - odel).abs();
619            if i > 1 && dd >= odd && ddel >= oddel {
620                break;
621            }
622            odd = dd;
623            oddel = ddel;
624        }
625        od = d;
626        odel = del;
627        i += 1;
628    }
629
630    if i >= IMAX {
631        warn += 4;
632    }
633
634    let ut = eraSxp_safe(d, &ust)?;
635    let ur = eraSxp_safe(ERFA_DC * (d * betsr + del), &pu)?;
636    let vfin = eraPpp_safe(&ur, &ut)?;
637    pv[1] = vfin;
638
639    Ok((pv, warn))
640}
641
642//----------------------------------------------------------------------
643// G30/sxp.c
644//----------------------------------------------------------------------
645// Scalar × p-vector.
646pub fn eraSxp_safe(s: f64, p: &[f64; 3]) -> ErfaResult<[f64; 3]> {
647    Ok([s * p[0], s * p[1], s * p[2]])
648}
649
650//----------------------------------------------------------------------
651// G30/sxpv.c
652//----------------------------------------------------------------------
653// Scalar × pv-vector.
654pub fn eraSxpv_safe(s: f64, pv: &[[f64; 3]; 2]) -> ErfaResult<[[f64; 3]; 2]> {
655    eraS2xpv_safe(s, s, pv)
656}
erfa_rust/G30_safe.rs

erfa_rust/
G30_safe.rs
