erfa_rust/
G25_safe.rs

1// G25
2//   pmat00.c → eraPmat00_safe
3//   pmat06.c → eraPmat06_safe
4//   pmat76.c → eraPmat76_safe
5//   pmp.c    → eraPmp_safe
6//   pmpx.c   → eraPmpx_safe
7//   pmsafe.c → eraPmsafe_safe
8//   pn.c     → eraPn_safe
9//   pn00.c   → eraPn00_safe
10//   pn00a.c  → eraPn00a_safe
11//   pn00b.c  → eraPn00b_safe
12
13use crate::G16_safe::eraFw2m_safe;
14use crate::G19_safe::eraIr_safe;
15use crate::G21_safe::eraNumat_safe;
16use crate::G22_safe::eraNut00a_safe;
17use crate::G23_safe::{eraNut00b_safe, eraObl80_safe};
18use crate::G24_safe::{eraPdp_safe, eraPfw06_safe, eraPm_safe};
19use crate::G26_safe::eraPr00_safe;
20use crate::G27_safe::eraPrec76_safe;
21use crate::G28_safe::{eraRxr_safe, eraRy_safe, eraRz_safe};
22use crate::G30_safe::{eraSeps_safe, eraStarpm_safe, eraSxp_safe};
23use crate::G35_safe::eraZp_safe;
24use crate::G6_safe::eraBp00_safe;
25use crate::H1_safe::{ERFA_AULT, ERFA_DAS2R, ERFA_DAU, ERFA_DAYSEC, ERFA_DJ00, ERFA_DJM, ERFA_DJY};
26
27pub type ErfaResult<T> = Result<T, ()>;
28
29// Precession matrix, IAU 2000 bias-precession model.
30pub fn eraPmat00_safe(date1: f64, date2: f64) -> ErfaResult<[[f64; 3]; 3]> {
31    let (_rb, _rp, rbp) = eraBp00_safe(date1, date2)?;
32    Ok(rbp)
33}
34
35// Precession matrix, IAU 2006 bias-precession model.
36pub fn eraPmat06_safe(date1: f64, date2: f64) -> ErfaResult<[[f64; 3]; 3]> {
37    let (gamb, phib, psib, epsa) = eraPfw06_safe(date1, date2)?;
38    let rbp = eraFw2m_safe(gamb, phib, psib, epsa)?;
39    Ok(rbp)
40}
41
42// Precession matrix, IAU 1976 model.
43pub fn eraPmat76_safe(date1: f64, date2: f64) -> ErfaResult<[[f64; 3]; 3]> {
44    let (zeta, z, theta) = eraPrec76_safe(ERFA_DJ00, 0.0, date1, date2)?;
45    let mut r = [[0.0_f64; 3]; 3];
46    eraIr_safe(&mut r)?;
47    eraRz_safe(-zeta, &mut r)?;
48    eraRy_safe(theta, &mut r)?;
49    eraRz_safe(-z, &mut r)?;
50    Ok(r)
51}
52
53// P-vector subtraction (a − b).
54pub fn eraPmp_safe(a: &[f64; 3], b: &[f64; 3]) -> ErfaResult<[f64; 3]> {
55    Ok([a[0] - b[0], a[1] - b[1], a[2] - b[2]])
56}
57
58// Proper motion and parallax: catalog position → updated unit vector.
59pub fn eraPmpx_safe(
60    rc: f64,
61    dc: f64,
62    pr: f64,
63    pd: f64,
64    px: f64,
65    rv: f64,
66    pmt: f64,
67    pob: &[f64; 3],
68) -> ErfaResult<[f64; 3]> {
69    // Unit conversions and constants.
70    const VF: f64 = ERFA_DAYSEC * ERFA_DJM / ERFA_DAU; // km/s → au/yr
71    const AULTY: f64 = ERFA_AULT / ERFA_DAYSEC / ERFA_DJY; // 1 au light-time (yr)
72
73    // Spherical to Cartesian for catalog direction.
74    let (sr, cr) = rc.sin_cos();
75    let (sd, cd) = dc.sin_cos();
76    let mut p = [cr * cd, sr * cd, sd];
77
78    // Effective interval with Roemer term.
79    let dt = pmt + eraPdp_safe(&p, pob)? * AULTY;
80
81    // Space motion (rad/yr) and parallax.
82    let pxr = px * ERFA_DAS2R;
83    let w = VF * rv * pxr;
84    let pdz = pd * p[2];
85    let pm = [
86        -pr * p[1] - pdz * cr + w * p[0],
87        pr * p[0] - pdz * sr + w * p[1],
88        pd * cd + w * p[2],
89    ];
90
91    // Apply motion and parallax.
92    for i in 0..3 {
93        p[i] += dt * pm[i] - pxr * pob[i];
94    }
95
96    // Normalize to unit vector.
97    let (_r, u) = eraPn_safe(&p)?;
98    Ok(u)
99}
100
101// Proper-motion propagation with safe distance override.
102pub fn eraPmsafe_safe(
103    ra1: f64,
104    dec1: f64,
105    pmr1: f64,
106    pmd1: f64,
107    px1: f64,
108    rv1: f64,
109    ep1a: f64,
110    ep1b: f64,
111    ep2a: f64,
112    ep2b: f64,
113) -> ErfaResult<((f64, f64, f64, f64, f64, f64), i32)> {
114    const PXMIN: f64 = 5.0e-7; // arcsec
115    const F: f64 = 326.0; // scale giving ~1% c max transverse speed
116
117    // Proper motion magnitude over one year (radians).
118    let pm = eraSeps_safe(ra1, dec1, ra1 + pmr1, dec1 + pmd1)?;
119
120    // Override small parallax if implausible given PM.
121    let mut jpx = 0;
122    let mut px1a = px1;
123    let pm_scaled = pm * F;
124    if px1a < pm_scaled {
125        px1a = pm_scaled;
126        jpx = 1;
127    }
128    if px1a < PXMIN {
129        px1a = PXMIN;
130        jpx = 1;
131    }
132
133    // Propagate star parameters.
134    let ((ra2, dec2, pmr2, pmd2, px2, rv2), mut j) =
135        eraStarpm_safe(ra1, dec1, pmr1, pmd1, px1a, rv1, ep1a, ep1b, ep2a, ep2b)?;
136
137    // Revise status to reflect parallax override if no error bits set.
138    if (j & 1) == 0 {
139        j += jpx;
140    }
141
142    Ok(((ra2, dec2, pmr2, pmd2, px2, rv2), j))
143}
144
145// Decompose p-vector into modulus and unit vector.
146pub fn eraPn_safe(p: &[f64; 3]) -> ErfaResult<(f64, [f64; 3])> {
147    let w = eraPm_safe(p)?;
148    if w == 0.0 {
149        let u = eraZp_safe();
150        Ok((w, u))
151    } else {
152        let u = eraSxp_safe(1.0 / w, p)?;
153        Ok((w, u))
154    }
155}
156
157// Precession-nutation, IAU 2000, with supplied dpsi,deps.
158pub fn eraPn00_safe(
159    date1: f64,
160    date2: f64,
161    dpsi: f64,
162    deps: f64,
163) -> ErfaResult<(
164    f64,           // epsa
165    [[f64; 3]; 3], // rb
166    [[f64; 3]; 3], // rp
167    [[f64; 3]; 3], // rbp
168    [[f64; 3]; 3], // rn
169    [[f64; 3]; 3], // rbpn
170)> {
171    // IAU 2000 precession-rate corrections.
172    let (_dpsipr, depspr) = eraPr00_safe(date1, date2)?;
173
174    // Mean obliquity corrected by depspr only.
175    let epsa = eraObl80_safe(date1, date2)? + depspr;
176
177    // Bias and precession matrices.
178    let (rb, rp, rbpw) = eraBp00_safe(date1, date2)?;
179    let rbp = rbpw;
180
181    // Nutation matrix built from given dpsi,deps (do not add dpsipr to dpsi).
182    let rnw = eraNumat_safe(epsa, dpsi, deps)?;
183    let rn = rnw;
184
185    // Bias-precession-nutation.
186    let rbpn = eraRxr_safe(&rnw, &rbpw)?;
187
188    Ok((epsa, rb, rp, rbp, rn, rbpn))
189}
190
191// Precession-nutation, IAU 2000A model.
192pub fn eraPn00a_safe(
193    date1: f64,
194    date2: f64,
195) -> ErfaResult<(
196    f64,           // dpsi
197    f64,           // deps
198    f64,           // epsa
199    [[f64; 3]; 3], // rb
200    [[f64; 3]; 3], // rp
201    [[f64; 3]; 3], // rbp
202    [[f64; 3]; 3], // rn
203    [[f64; 3]; 3], // rbpn
204)> {
205    let (dpsi, deps) = eraNut00a_safe(date1, date2)?;
206    let (epsa, rb, rp, rbp, rn, rbpn) = eraPn00_safe(date1, date2, dpsi, deps)?;
207    Ok((dpsi, deps, epsa, rb, rp, rbp, rn, rbpn))
208}
209
210// Precession-nutation, IAU 2000B model.
211pub fn eraPn00b_safe(
212    date1: f64,
213    date2: f64,
214) -> ErfaResult<(
215    f64,           // dpsi
216    f64,           // deps
217    f64,           // epsa
218    [[f64; 3]; 3], // rb
219    [[f64; 3]; 3], // rp
220    [[f64; 3]; 3], // rbp
221    [[f64; 3]; 3], // rn
222    [[f64; 3]; 3], // rbpn
223)> {
224    let (dpsi, deps) = eraNut00b_safe(date1, date2)?;
225    let (epsa, rb, rp, rbp, rn, rbpn) = eraPn00_safe(date1, date2, dpsi, deps)?;
226    Ok((dpsi, deps, epsa, rb, rp, rbp, rn, rbpn))
227}
erfa_rust/G25_safe.rs

erfa_rust/
G25_safe.rs
