erfa_rust/
G6_safe.rs

1// G6
2//   bi00.c   → eraBi00_safe
3//   bp00.c   → eraBp00_safe
4//   bp06.c   → eraBp06_safe
5//   bpn2xy.c → eraBpn2xy_safe
6
7use crate::H1_safe::{ERFA_DAS2R, ERFA_DJ00, ERFA_DJC, ERFA_DJM0, ERFA_DJM00};
8
9use crate::G16_safe::eraFw2m_safe;
10use crate::G19_safe::eraIr_safe;
11use crate::G24_safe::eraPfw06_safe;
12use crate::G25_safe::eraPmat06_safe;
13use crate::G26_safe::eraPr00_safe;
14use crate::G28_safe::{eraRx_safe, eraRxr_safe, eraRy_safe, eraRz_safe};
15use crate::G33_safe::eraTr_safe;
16use crate::G8_safe::eraCr_safe;
17
18pub type ErfaResult<T> = Result<T, ()>;
19
20// eraBi00_safe → bi00.c
21// Frame-bias corrections (ICRS → J2000.0), radians; returns (dpsibi, depsbi, dra).
22pub fn eraBi00_safe() -> ErfaResult<(f64, f64, f64)> {
23    const DPBIAS: f64 = -0.041_775 * ERFA_DAS2R;
24    const DEBIAS: f64 = -0.006_819_2 * ERFA_DAS2R;
25    const DRA0: f64 = -0.0146 * ERFA_DAS2R;
26    Ok((DPBIAS, DEBIAS, DRA0))
27}
28
29// eraBp00_safe → bp00.c
30// Frame bias (B), precession (P) and product (BP), IAU 2000; returns (rb, rp, rbp).
31pub fn eraBp00_safe(
32    date1: f64,
33    date2: f64,
34) -> ErfaResult<([[f64; 3]; 3], [[f64; 3]; 3], [[f64; 3]; 3])> {
35    const EPS0: f64 = 84_381.448 * ERFA_DAS2R;
36    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
37
38    let (dpsibi, depsbi, dra0) = eraBi00_safe()?;
39
40    let psia77 = (5038.7784 + (-1.07259 + -0.001_147 * t) * t) * t * ERFA_DAS2R;
41    let oma77 = EPS0 + (0.05127 + -0.007_726 * t) * t * t * ERFA_DAS2R;
42    let chia = (10.5526 + (-2.38064 + -0.001_125 * t) * t) * t * ERFA_DAS2R;
43
44    let (dpsipr, depspr) = eraPr00_safe(date1, date2)?;
45    let psia = psia77 + dpsipr;
46    let oma = oma77 + depspr;
47
48    let mut rbw = [[0.0_f64; 3]; 3];
49    let mut rb = [[0.0_f64; 3]; 3];
50    let mut rp = [[0.0_f64; 3]; 3];
51
52    // Bias: GCRS → J2000.0
53    eraIr_safe(&mut rbw)?;
54    eraRz_safe(dra0, &mut rbw)?;
55    eraRy_safe(dpsibi * EPS0.sin(), &mut rbw)?;
56    eraRx_safe(-depsbi, &mut rbw)?;
57    eraCr_safe(&rbw, &mut rb)?;
58
59    // Precession: J2000.0 → mean of date
60    eraIr_safe(&mut rp)?;
61    eraRx_safe(EPS0, &mut rp)?;
62    eraRz_safe(-psia, &mut rp)?;
63    eraRx_safe(-oma, &mut rp)?;
64    eraRz_safe(chia, &mut rp)?;
65
66    let rbp = eraRxr_safe(&rp, &rbw)?;
67    Ok((rb, rp, rbp))
68}
69
70// eraBp06_safe → bp06.c
71// Frame bias (B), precession (P) and product (BP), IAU 2006; returns (rb, rp, rbp).
72pub fn eraBp06_safe(
73    date1: f64,
74    date2: f64,
75) -> ErfaResult<([[f64; 3]; 3], [[f64; 3]; 3], [[f64; 3]; 3])> {
76    // Bias matrix from FukushimaWilliams angles at J2000.0
77    let (gamb, phib, psib, epsa) = eraPfw06_safe(ERFA_DJM0, ERFA_DJM00)?;
78    let rb = eraFw2m_safe(gamb, phib, psib, epsa)?;
79
80    // Precession matrix for given date (IAU 2006)
81    let rbpw = eraPmat06_safe(date1, date2)?;
82
83    // rp = P × Tr(B)
84    let rbt = eraTr_safe(&rb)?;
85    let rp = eraRxr_safe(&rbpw, &rbt)?;
86
87    // rbp = P
88    let rbp = rbpw;
89
90    Ok((rb, rp, rbp))
91}
92
93// eraBpn2xy_safe → bpn2xy.c
94// Extract CIP X,Y from rbpn matrix; returns (x, y) from row 2.
95pub fn eraBpn2xy_safe(rbpn: &[[f64; 3]; 3]) -> ErfaResult<(f64, f64)> {
96    let x = rbpn[2][0];
97    let y = rbpn[2][1];
98    Ok((x, y))
99}