erfa_rust/
G2_safe.rs

1// G2
2//   apco.c    → eraApco_safe
3//   apco13.c  → eraApco13_safe
4//   apcs.c    → eraApcs_safe
5//   apcs13.c  → eraApcs13_safe
6//   aper.c    → eraAper_safe
7
8
9use crate::G11_safe::eraEors_safe;
10use crate::G13_safe::eraEpv00_safe;
11use crate::G14_safe::eraEra00_safe;
12use crate::G19_safe::eraIr_safe;
13use crate::G1_safe::eraAnpm_safe;
14use crate::G25_safe::eraPn_safe;
15use crate::G26_safe::eraPnm06a_safe;
16use crate::G27_safe::eraPvtob_safe;
17use crate::G28_safe::{eraRefco_safe, eraRx_safe, eraRy_safe, eraRz_safe};
18use crate::G30_safe::{eraS06_safe, eraSp00_safe};
19use crate::G32_safe::eraTaitt_safe;
20use crate::G33_safe::{eraTrxpv_safe, eraUtctai_safe, eraUtcut1_safe};
21use crate::G6_safe::eraBpn2xy_safe;
22use crate::G7_safe::eraC2ixys_safe;
23use crate::G8_safe::{eraCp_safe, eraCr_safe};
24use crate::H1_safe::{eraASTROM, ERFA_AULT, ERFA_DAU, ERFA_DAYSEC, ERFA_DJ00, ERFA_DJY};
25
26pub type ErfaResult<T> = Result<T, ()>;
27
28// G2/apco.c
29// Prepare star-independent astrometry parameters for a terrestrial observer.
30pub fn eraApco_safe(
31    date1: f64,
32    date2: f64,
33    ebpv: &[[f64; 3]; 2],
34    ehp: &[f64; 3],
35    x: f64,
36    y: f64,
37    s: f64,
38    theta: f64,
39    elong: f64,
40    phi: f64,
41    hm: f64,
42    xp: f64,
43    yp: f64,
44    sp: f64,
45    refa: f64,
46    refb: f64,
47    astrom: &mut eraASTROM,
48) -> ErfaResult<()> {
49    // Build rotation matrix for local parameters.
50    let mut r = [[0.0_f64; 3]; 3];
51    eraIr_safe(&mut r)?;
52    eraRz_safe(theta + sp, &mut r)?;
53    eraRy_safe(-xp, &mut r)?;
54    eraRx_safe(-yp, &mut r)?;
55    eraRz_safe(elong, &mut r)?;
56
57    // Local Earth rotation angle.
58    let mut a = r[0][0];
59    let mut b = r[0][1];
60    let eral = if a != 0.0 || b != 0.0 {
61        b.atan2(a)
62    } else {
63        0.0
64    };
65    astrom.eral = eral;
66
67    // Polar motion with respect to local meridian.
68    a = r[0][0];
69    let c = r[0][2];
70    astrom.xpl = c.atan2((a * a + b * b).sqrt());
71    a = r[1][2];
72    b = r[2][2];
73    astrom.ypl = if a != 0.0 || b != 0.0 {
74        -a.atan2(b)
75    } else {
76        0.0
77    };
78
79    // Adjusted longitude.
80    astrom.along = eraAnpm_safe(eral - theta)?;
81
82    // Latitude functions.
83    astrom.sphi = phi.sin();
84    astrom.cphi = phi.cos();
85
86    // Refraction constants.
87    astrom.refa = refa;
88    astrom.refb = refb;
89
90    // Disable diurnal aberration step.
91    astrom.diurab = 0.0;
92
93    // CIO-based C2I matrix (overwrite r with returned matrix).
94    let r_c2i = eraC2ixys_safe(x, y, s)?;
95    r = r_c2i;
96
97    // Observer geocentric PV in CIRS, then rotate into GCRS.
98    let pvc = eraPvtob_safe(elong, phi, hm, xp, yp, sp, theta)?;
99    let pv = eraTrxpv_safe(&r, &pvc)?;
100
101    // ICRS <-> GCRS parameters.
102    eraApcs_safe(date1, date2, &pv, ebpv, ehp, astrom)?;
103
104    // Store the C2I (BPN) matrix.
105    eraCr_safe(&r, &mut astrom.bpn)?;
106
107    Ok(())
108}
109
110// G2/apco13.c
111// Prepare astrometry parameters from UTC; returns (eo, j) where j=0 or +1.
112pub fn eraApco13_safe(
113    utc1: f64,
114    utc2: f64,
115    dut1: f64,
116    elong: f64,
117    phi: f64,
118    hm: f64,
119    xp: f64,
120    yp: f64,
121    phpa: f64,
122    tc: f64,
123    rh: f64,
124    wl: f64,
125    astrom: &mut eraASTROM,
126) -> ErfaResult<(f64, i32)> {
127    // UTC→TAI and TAI→TT; UTC→UT1.
128    let ((tai1, tai2), j_utctai) = eraUtctai_safe(utc1, utc2)?;
129    if j_utctai < 0 {
130        return Err(());
131    }
132    let ((tt1, tt2), _j_tai_tt) = eraTaitt_safe(tai1, tai2)?;
133    let ((ut11, ut12), j_ut1) = eraUtcut1_safe(utc1, utc2, dut1)?;
134    if j_ut1 < 0 {
135        return Err(());
136    }
137
138    // Earth ephemeris, CIP/CIO, refraction.
139    let (ehpv, ebpv, _jstat) = eraEpv00_safe(tt1, tt2)?;
140    let r = eraPnm06a_safe(tt1, tt2)?;
141    let (x, y) = eraBpn2xy_safe(&r)?;
142    let s = eraS06_safe(tt1, tt2, x, y)?;
143    let theta = eraEra00_safe(ut11, ut12)?;
144    let sp = eraSp00_safe(tt1, tt2)?;
145    let (refa, refb) = eraRefco_safe(phpa, tc, rh, wl)?;
146
147    // Assemble star-independent astrometry parameters.
148    eraApco_safe(
149        tt1, tt2, &ebpv, &ehpv[0], x, y, s, theta, elong, phi, hm, xp, yp, sp, refa, refb, astrom,
150    )?;
151
152    // Equation of the origins.
153    let eo = eraEors_safe(&r, s)?;
154
155    // Return EO and UT1 status (0 or +1).
156    Ok((eo, j_ut1))
157}
158
159// G2/apcs.c
160// Compute ICRS↔GCRS parameters for a space observer.
161pub fn eraApcs_safe(
162    date1: f64,
163    date2: f64,
164    pv: &[[f64; 3]; 2],   // (m, m/s)
165    ebpv: &[[f64; 3]; 2], // (au, au/d)
166    ehp: &[f64; 3],       // (au)
167    astrom: &mut eraASTROM,
168) -> ErfaResult<()> {
169    // Unit conversions.
170    const AUDMS: f64 = ERFA_DAU / ERFA_DAYSEC; // au/day to m/s
171    const CR: f64 = ERFA_AULT / ERFA_DAYSEC; // speed of light in au/day
172
173    // Time since reference epoch, years (for proper motion).
174    astrom.pmt = ((date1 - ERFA_DJ00) + date2) / ERFA_DJY;
175
176    // Adjust Earth ephemeris to observer.
177    let mut pb = [0.0_f64; 3];
178    let mut vb = [0.0_f64; 3];
179    let mut ph = [0.0_f64; 3];
180    for i in 0..3 {
181        let dp = pv[0][i] / ERFA_DAU;
182        let dv = pv[1][i] / AUDMS;
183        pb[i] = ebpv[0][i] + dp;
184        vb[i] = ebpv[1][i] + dv;
185        ph[i] = ehp[i] + dp;
186    }
187
188    // Barycentric position of observer (au).
189    eraCp_safe(&pb, &mut astrom.eb)?;
190
191    // Heliocentric direction (unit) and distance (au).
192    let (em, eh) = eraPn_safe(&ph)?;
193    astrom.em = em;
194    astrom.eh = eh;
195
196    // Barycentric velocity in units of c and reciprocal Lorentz factor.
197    let mut v2 = 0.0;
198    for i in 0..3 {
199        let w = vb[i] * CR;
200        astrom.v[i] = w;
201        v2 += w * w;
202    }
203    astrom.bm1 = (1.0 - v2).sqrt();
204
205    // Reset the NPB matrix.
206    eraIr_safe(&mut astrom.bpn)?;
207
208    Ok(())
209}
210
211// G2/apcs13.c
212// As eraApcs_safe but Earth ephemerides computed internally.
213pub fn eraApcs13_safe(
214    date1: f64,
215    date2: f64,
216    pv: &[[f64; 3]; 2], // (m, m/s)
217    astrom: &mut eraASTROM,
218) -> ErfaResult<()> {
219    // Earth barycentric & heliocentric PV (au, au/day).
220    let (ehpv, ebpv, _jstat) = eraEpv00_safe(date1, date2)?;
221
222    // Compute the star-independent astrometry parameters.
223    eraApcs_safe(date1, date2, pv, &ebpv, &ehpv[0], astrom)
224}
225
226// G2/aper.c
227// Update local apparent sidereal angle.
228pub fn eraAper_safe(theta: f64, astrom: &mut eraASTROM) -> ErfaResult<()> {
229    astrom.eral = theta + astrom.along;
230    Ok(())
231}
erfa_rust/G2_safe.rs

erfa_rust/
G2_safe.rs
