erfa_rust/
G18_safe.rs

1// G18
2//   h2fk5.c  → eraH2fk5_safe
3//   hd2ae.c  → eraHd2ae_safe
4//   hd2pa.c  → eraHd2pa_safe
5//   hfk5z.c  → eraHfk5z_safe
6
7use crate::G15_safe::eraFk5hip_safe;
8use crate::G1_safe::eraAnp_safe;
9use crate::G25_safe::eraPmp_safe;
10use crate::G27_safe::{eraPv2s_safe, eraPvstar_safe, eraPxp_safe};
11use crate::G28_safe::{eraRv2m_safe, eraRxp_safe, eraRxr_safe};
12use crate::G29_safe::eraS2c_safe;
13use crate::G30_safe::{eraStarpv_safe, eraSxp_safe};
14use crate::G33_safe::eraTrxp_safe;
15use crate::H1_safe::{ERFA_D2PI, ERFA_DJ00, ERFA_DJY};
16
17pub type ErfaResult<T> = Result<T, ()>;
18
19//----------------------------------------------------------------------
20// G18/h2fk5.c → eraH2fk5_safe
21//----------------------------------------------------------------------
22// Convert Hipparcos catalog parameters to FK5 J2000 parameters.
23pub fn eraH2fk5_safe(
24    rh: f64,
25    dh: f64,
26    drh: f64,
27    ddh: f64,
28    pxh: f64,
29    rvh: f64,
30) -> ErfaResult<(f64, f64, f64, f64, f64, f64)> {
31    // Hipparcos barycentric position/velocity (normalized)
32    let (pvh, _warn) = eraStarpv_safe(rh, dh, drh, ddh, pxh, rvh)?;
33
34    // FK5↔Hipparcos rotation and spin (spin in rad/yr)
35    let (r5h_m, mut s5h) = eraFk5hip_safe()?;
36
37    // Spin: yr⁻¹ → day⁻¹
38    for v in &mut s5h {
39        *v /= 365.25;
40    }
41
42    // Spin expressed in Hipparcos frame: sh = R * s5h
43    let sh = eraRxp_safe(&r5h_m, &s5h)?;
44
45    // De-orient position: pv5.pos = R^T * pvh.pos
46    let mut pv5 = [[0.0_f64; 3]; 2];
47    pv5[0] = eraTrxp_safe(&r5h_m, &pvh[0])?;
48
49    // Extra motion from spin: wxp = pvh.pos × sh
50    let wxp = eraPxp_safe(&pvh[0], &sh)?;
51
52    // Remove spin component from velocity: vv = pvh.vel − wxp
53    let vv = eraPmp_safe(&pvh[1], &wxp)?;
54
55    // De-orient velocity: pv5.vel = R^T * vv
56    pv5[1] = eraTrxp_safe(&r5h_m, &vv)?;
57
58    // FK5 pv-vector → catalog fields
59    let ((r5, d5, dr5, dd5, px5, rv5), _j) = eraPvstar_safe(&pv5)?;
60    Ok((r5, d5, dr5, dd5, px5, rv5))
61}
62
63//----------------------------------------------------------------------
64// G18/hd2ae.c → eraHd2ae_safe
65//----------------------------------------------------------------------
66// Convert hour angle/declination and latitude to azimuth/elevation.
67pub fn eraHd2ae_safe(ha: f64, dec: f64, phi: f64) -> ErfaResult<(f64, f64)> {
68    let (sh, ch) = ha.sin_cos();
69    let (sd, cd) = dec.sin_cos();
70    let (sp, cp) = phi.sin_cos();
71
72    // Horizon-system unit vector
73    let x = -ch * cd * sp + sd * cp;
74    let y = -sh * cd;
75    let z = ch * cd * cp + sd * sp;
76
77    // Spherical coordinates
78    let r = (x * x + y * y).sqrt();
79    let mut a = if r != 0.0 { y.atan2(x) } else { 0.0 };
80    if a < 0.0 {
81        a += ERFA_D2PI;
82    }
83
84    let az = a;
85    let el = z.atan2(r);
86    Ok((az, el))
87}
88
89//----------------------------------------------------------------------
90// G18/hd2pa.c → eraHd2pa_safe
91//----------------------------------------------------------------------
92// Compute parallactic angle for hour angle/declination and latitude.
93pub fn eraHd2pa_safe(ha: f64, dec: f64, phi: f64) -> ErfaResult<f64> {
94    let cp = phi.cos();
95    let sqsz = cp * ha.sin();
96    let cqsz = phi.sin() * dec.cos() - cp * dec.sin() * ha.cos();
97    let pa = if sqsz != 0.0 || cqsz != 0.0 {
98        sqsz.atan2(cqsz)
99    } else {
100        0.0
101    };
102    Ok(pa)
103}
104
105//----------------------------------------------------------------------
106// G18/hfk5z.c → eraHfk5z_safe
107//----------------------------------------------------------------------
108// Convert Hipparcos position (zero proper motion) to FK5 at given date.
109pub fn eraHfk5z_safe(rh: f64, dh: f64, date1: f64, date2: f64) -> ErfaResult<(f64, f64, f64, f64)> {
110    // Interval from J2000.0 in Julian years
111    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJY;
112
113    // Hipparcos unit vector
114    let ph = eraS2c_safe(rh, dh)?;
115
116    // FK5↔Hipparcos rotation and spin (spin in rad/yr)
117    let (r5h_m, s5h) = eraFk5hip_safe()?;
118
119    // Spin rotated into Hipparcos: sh = R * s5h
120    let sh = eraRxp_safe(&r5h_m, &s5h)?;
121
122    // Accumulated spin over interval (vector): vst = t  s5h
123    let vst = eraSxp_safe(t, &s5h)?;
124
125    // Vector → rotation matrix
126    let rst = eraRv2m_safe(&vst)?;
127
128    // Total rotation: r5ht = R × rst
129    let r5ht = eraRxr_safe(&r5h_m, &rst)?;
130
131    // De-orient + de-spin position: pv5e.pos = r5ht^T * ph
132    let mut pv5e = [[0.0_f64; 3]; 2];
133    pv5e[0] = eraTrxp_safe(&r5ht, &ph)?;
134
135    // Space motion due to spin: vv = sh × ph
136    let vv = eraPxp_safe(&sh, &ph)?;
137
138    // De-orient + de-spin velocity: pv5e.vel = r5ht^T * vv
139    pv5e[1] = eraTrxp_safe(&r5ht, &vv)?;
140
141    // PV → spherical (theta, phi, r, td, pd, rd)
142    let (w, d5, _r, dr5, dd5, _rd) = eraPv2s_safe(&pv5e)?;
143    let r5 = eraAnp_safe(w)?;
144    Ok((r5, d5, dr5, dd5))
145}
erfa_rust/G18_safe.rs

erfa_rust/
G18_safe.rs
