erfa_rust/
G27_safe.rs

1// G27
2//   prec76.c  → eraPrec76_safe
3//   pv2p.c    → eraPv2p_safe
4//   pv2s.c    → eraPv2s_safe
5//   pvdpv.c   → eraPvdpv_safe
6//   pvm.c     → eraPvm_safe
7//   pvmpv.c   → eraPvmpv_safe
8//   pvppv.c   → eraPvppv_safe
9//   pvstar.c  → eraPvstar_safe
10//   pvtob.c   → eraPvtob_safe
11//   pvu.c     → eraPvu_safe
12//   pvup.c    → eraPvup_safe
13//   pvxpv.c   → eraPvxpv_safe
14//   pxp.c     → eraPxp_safe
15use crate::H1_safe::{
16    ERFA_D2PI, ERFA_DAS2R, ERFA_DAU, ERFA_DAYSEC, ERFA_DC, ERFA_DJ00, ERFA_DJC, ERFA_DJY,
17    ERFA_DR2AS,
18};
19
20use crate::G17_safe::eraGd2gc_safe;
21use crate::G1_safe::eraAnp_safe;
22use crate::G24_safe::{eraPdp_safe, eraPm_safe};
23use crate::G25_safe::{eraPmp_safe, eraPn_safe};
24use crate::G26_safe::{eraPom00_safe, eraPpp_safe, eraPpsp_safe};
25use crate::G30_safe::eraSxp_safe;
26use crate::G33_safe::eraTrxp_safe;
27use crate::G8_safe::{eraCp_safe, eraCpv_safe};
28
29pub type ErfaResult<T> = Result<T, ()>;
30
31// Compute IAU 1976 precession angles (zeta, z, theta) between two epochs.
32pub fn eraPrec76_safe(
33    date01: f64,
34    date02: f64,
35    date11: f64,
36    date12: f64,
37) -> ErfaResult<(f64, f64, f64)> {
38    let t0 = ((date01 - ERFA_DJ00) + date02) / ERFA_DJC;
39    let t = ((date11 - date01) + (date12 - date02)) / ERFA_DJC;
40    let tas2r = t * ERFA_DAS2R;
41    let w = 2306.2181 + (1.39656 - 0.000139 * t0) * t0;
42    let zeta = (w + ((0.30188 - 0.000344 * t0) + 0.017998 * t) * t) * tas2r;
43    let z = (w + ((1.09468 + 0.000066 * t0) + 0.018203 * t) * t) * tas2r;
44    let theta = ((2004.3109 + (-0.85330 - 0.000217 * t0) * t0)
45        + ((-0.42665 - 0.000217 * t0) - 0.041833 * t) * t)
46        * tas2r;
47    Ok((zeta, z, theta))
48}
49
50// Extract position component from pv-vector.
51pub fn eraPv2p_safe(pv: &[[f64; 3]; 2]) -> ErfaResult<[f64; 3]> {
52    Ok(pv[0])
53}
54
55// Convert pv-vector to spherical angles and their rates.
56pub fn eraPv2s_safe(pv: &[[f64; 3]; 2]) -> ErfaResult<(f64, f64, f64, f64, f64, f64)> {
57    let mut x = pv[0][0];
58    let mut y = pv[0][1];
59    let mut z = pv[0][2];
60    let xd = pv[1][0];
61    let yd = pv[1][1];
62    let zd = pv[1][2];
63
64    let mut rxy2 = x * x + y * y;
65    let mut r2 = rxy2 + z * z;
66    let rtrue = r2.sqrt();
67
68    let rw = if rtrue == 0.0 {
69        x = xd;
70        y = yd;
71        z = zd;
72        rxy2 = x * x + y * y;
73        r2 = rxy2 + z * z;
74        r2.sqrt()
75    } else {
76        rtrue
77    };
78
79    let rxy = rxy2.sqrt();
80    let xyp = x * xd + y * yd;
81    let (theta, phi, td, pd) = if rxy2 != 0.0 {
82        let theta = y.atan2(x);
83        let phi = z.atan2(rxy);
84        let td = (x * yd - y * xd) / rxy2;
85        let pd = (zd * rxy2 - z * xyp) / (r2 * rxy);
86        (theta, phi, td, pd)
87    } else {
88        let theta = 0.0;
89        let phi = if z != 0.0 { z.atan2(rxy) } else { 0.0 };
90        let td = 0.0;
91        let pd = 0.0;
92        (theta, phi, td, pd)
93    };
94    let r = rtrue;
95    let rd = if rw != 0.0 { (xyp + z * zd) / rw } else { 0.0 };
96
97    Ok((theta, phi, r, td, pd, rd))
98}
99
100// Dot product of two pv-vectors and its time derivative.
101pub fn eraPvdpv_safe(a: &[[f64; 3]; 2], b: &[[f64; 3]; 2]) -> ErfaResult<(f64, f64)> {
102    let p_dot = eraPdp_safe(&a[0], &b[0])?;
103    let adbd = eraPdp_safe(&a[0], &b[1])?;
104    let addb = eraPdp_safe(&a[1], &b[0])?;
105    let d_dot = adbd + addb;
106    Ok((p_dot, d_dot))
107}
108
109// Moduli of position and velocity from a pv-vector.
110pub fn eraPvm_safe(pv: &[[f64; 3]; 2]) -> ErfaResult<(f64, f64)> {
111    let r = eraPm_safe(&pv[0])?;
112    let s = eraPm_safe(&pv[1])?;
113    Ok((r, s))
114}
115
116// Subtract two pv-vectors component-wise.
117pub fn eraPvmpv_safe(a: &[[f64; 3]; 2], b: &[[f64; 3]; 2]) -> ErfaResult<[[f64; 3]; 2]> {
118    let p = eraPmp_safe(&a[0], &b[0])?;
119    let v = eraPmp_safe(&a[1], &b[1])?;
120    Ok([p, v])
121}
122
123// Add two pv-vectors component-wise.
124pub fn eraPvppv_safe(a: &[[f64; 3]; 2], b: &[[f64; 3]; 2]) -> ErfaResult<[[f64; 3]; 2]> {
125    let p = eraPpp_safe(&a[0], &b[0])?;
126    let v = eraPpp_safe(&a[1], &b[1])?;
127    Ok([p, v])
128}
129
130// Convert pv-vector to catalog parameters (ra,dec,pmr,pmd,px,rv); j=-1 superluminal, -2 null pos.
131pub fn eraPvstar_safe(pv: &[[f64; 3]; 2]) -> ErfaResult<((f64, f64, f64, f64, f64, f64), i32)> {
132    let (_r, pu) = eraPn_safe(&pv[0])?;
133    let vr = eraPdp_safe(&pu, &pv[1])?;
134    let ur = eraSxp_safe(vr, &pu)?;
135    let ut = eraPmp_safe(&pv[1], &ur)?;
136    let vt = eraPm_safe(&ut)?;
137
138    let bett = vt / ERFA_DC;
139    let betr = vr / ERFA_DC;
140    let d = 1.0 + betr;
141    let w = betr * betr + bett * bett;
142    if d == 0.0 || w > 1.0 {
143        return Ok(((0.0, 0.0, 0.0, 0.0, 0.0, 0.0), -1));
144    }
145    let del = -w / ((1.0 - w).sqrt() + 1.0);
146
147    let ust = eraSxp_safe(1.0 / d, &ut)?;
148    let usr = eraSxp_safe(ERFA_DC * (betr - del) / d, &pu)?;
149    let vnew = eraPpp_safe(&usr, &ust)?;
150    let pv_mod = [pv[0], vnew];
151
152    let (a, dec, r_out, rad, decd, rd) = eraPv2s_safe(&pv_mod)?;
153    if r_out == 0.0 {
154        return Ok(((0.0, 0.0, 0.0, 0.0, 0.0, 0.0), -2));
155    }
156
157    let ra = eraAnp_safe(a)?;
158    let pmr = rad * ERFA_DJY;
159    let pmd = decd * ERFA_DJY;
160    let px = ERFA_DR2AS / r_out;
161    let rv = 1e-3 * rd * ERFA_DAU / ERFA_DAYSEC;
162
163    Ok(((ra, dec, pmr, pmd, px, rv), 0))
164}
165
166// Observer geocentric position/velocity from site geodetic coordinates.
167pub fn eraPvtob_safe(
168    elong: f64,
169    phi: f64,
170    hm: f64,
171    xp: f64,
172    yp: f64,
173    sp: f64,
174    theta: f64,
175) -> ErfaResult<[[f64; 3]; 2]> {
176    const OM: f64 = 1.002_737_811_911_354_48 * ERFA_D2PI / ERFA_DAYSEC;
177
178    let (xyzm, _j) = eraGd2gc_safe(1, elong, phi, hm)?;
179    let rpm = eraPom00_safe(xp, yp, sp)?;
180    let xyz = eraTrxp_safe(&rpm, &xyzm)?;
181
182    let x = xyz[0];
183    let y = xyz[1];
184    let z = xyz[2];
185    let (s, c) = theta.sin_cos();
186
187    let px = c * x - s * y;
188    let py = s * x + c * y;
189    let pz = z;
190
191    let vx = OM * (-s * x - c * y);
192    let vy = OM * (c * x - s * y);
193    let vz = 0.0;
194
195    Ok([[px, py, pz], [vx, vy, vz]])
196}
197
198// Advance a pv-vector by time dt (same units as velocity denominator).
199pub fn eraPvu_safe(dt: f64, pv: &[[f64; 3]; 2]) -> ErfaResult<[[f64; 3]; 2]> {
200    let p = eraPpsp_safe(&pv[0], dt, &pv[1])?;
201    let mut v = [0.0_f64; 3];
202    eraCp_safe(&pv[1], &mut v)?;
203    Ok([p, v])
204}
205
206// Advance position only by time dt using a pv-vector.
207pub fn eraPvup_safe(dt: f64, pv: &[[f64; 3]; 2]) -> ErfaResult<[f64; 3]> {
208    Ok([
209        pv[0][0] + dt * pv[1][0],
210        pv[0][1] + dt * pv[1][1],
211        pv[0][2] + dt * pv[1][2],
212    ])
213}
214
215// Cross product of two pv-vectors, returning position and its derivative.
216pub fn eraPvxpv_safe(a: &[[f64; 3]; 2], b: &[[f64; 3]; 2]) -> ErfaResult<[[f64; 3]; 2]> {
217    let mut wa = [[0.0_f64; 3]; 2];
218    let mut wb = [[0.0_f64; 3]; 2];
219    eraCpv_safe(a, &mut wa)?;
220    eraCpv_safe(b, &mut wb)?;
221
222    let p = eraPxp_safe(&wa[0], &wb[0])?;
223    let axbd = eraPxp_safe(&wa[0], &wb[1])?;
224    let adxb = eraPxp_safe(&wa[1], &wb[0])?;
225    let v = eraPpp_safe(&axbd, &adxb)?;
226    Ok([p, v])
227}
228
229// 3D vector cross product.
230pub fn eraPxp_safe(a: &[f64; 3], b: &[f64; 3]) -> ErfaResult<[f64; 3]> {
231    let xa = a[0];
232    let ya = a[1];
233    let za = a[2];
234    let xb = b[0];
235    let yb = b[1];
236    let zb = b[2];
237
238    Ok([ya * zb - za * yb, za * xb - xa * zb, xa * yb - ya * xb])
239}
erfa_rust/G27_safe.rs

erfa_rust/
G27_safe.rs
