erfa_rust/
G33_safe.rs

1// G33
2//   tr.c      → eraTr_safe
3//   trxp.c    → eraTrxp_safe
4//   trxpv.c   → eraTrxpv_safe
5//   tttai.c   → eraTttai_safe
6//   tttcg.c   → eraTttcg_safe
7//   tttdb.c   → eraTttdb_safe
8//   ttut1.c   → eraTtut1_safe
9//   ut1tai.c  → eraUt1tai_safe
10//   ut1tt.c   → eraUt1tt_safe
11//   ut1utc.c  → eraUt1utc_safe
12//   utctai.c  → eraUtctai_safe
13//   utcut1.c  → eraUtcut1_safe
14
15use crate::G19_safe::eraJd2cal_safe;
16use crate::G28_safe::{eraRxp_safe, eraRxpv_safe};
17use crate::G32_safe::eraTaiut1_safe;
18use crate::G8_safe::{eraCal2jd_safe, eraCr_safe};
19use crate::G9_safe::eraDat_safe;
20use crate::H1_safe::{ERFA_DAYSEC, ERFA_DJM0, ERFA_DJM77, ERFA_ELG, ERFA_TTMTAI};
21
22pub type ErfaResult<T> = Result<T, ()>;
23
24// Transpose 3×3 matrix.
25pub fn eraTr_safe(r: &[[f64; 3]; 3]) -> ErfaResult<[[f64; 3]; 3]> {
26    let mut wm = [[0.0_f64; 3]; 3];
27    for i in 0..3 {
28        for j in 0..3 {
29            wm[i][j] = r[j][i];
30        }
31    }
32    // Copy to result (preserve original helper usage semantics)
33    let mut rt = [[0.0_f64; 3]; 3];
34    eraCr_safe(&wm, &mut rt)?;
35    Ok(rt)
36}
37
38// Transpose(r) × p-vector.
39pub fn eraTrxp_safe(r: &[[f64; 3]; 3], p: &[f64; 3]) -> ErfaResult<[f64; 3]> {
40    let tr = eraTr_safe(r)?;
41    eraRxp_safe(&tr, p)
42}
43
44// Transpose(r) × pv-vector.
45pub fn eraTrxpv_safe(r: &[[f64; 3]; 3], pv: &[[f64; 3]; 2]) -> ErfaResult<[[f64; 3]; 2]> {
46    let tr = eraTr_safe(r)?;
47    eraRxpv_safe(&tr, pv)
48}
49
50// TT → TAI.
51pub fn eraTttai_safe(tt1: f64, tt2: f64) -> ErfaResult<((f64, f64), i32)> {
52    let dtat = ERFA_TTMTAI / ERFA_DAYSEC;
53    let (tai1, tai2) = if tt1.abs() > tt2.abs() {
54        (tt1, tt2 - dtat)
55    } else {
56        (tt1 - dtat, tt2)
57    };
58    Ok(((tai1, tai2), 0))
59}
60
61// TT → TCG.
62pub fn eraTttcg_safe(tt1: f64, tt2: f64) -> ErfaResult<((f64, f64), i32)> {
63    const T77T: f64 = ERFA_DJM77 + ERFA_TTMTAI / ERFA_DAYSEC; // 1977-Jan-1 00:00:32.184 TT
64    const ELGG: f64 = ERFA_ELG / (1.0 - ERFA_ELG); // TT→TCG rate
65
66    let (tcg1, tcg2) = if tt1.abs() > tt2.abs() {
67        (tt1, tt2 + ((tt1 - ERFA_DJM0) + (tt2 - T77T)) * ELGG)
68    } else {
69        (tt1 + ((tt2 - ERFA_DJM0) + (tt1 - T77T)) * ELGG, tt2)
70    };
71    Ok(((tcg1, tcg2), 0))
72}
73
74// TT → TDB using caller-supplied dtr = TDB−TT seconds.
75pub fn eraTttdb_safe(tt1: f64, tt2: f64, dtr: f64) -> ErfaResult<((f64, f64), i32)> {
76    let dtrd = dtr / ERFA_DAYSEC;
77    let (tdb1, tdb2) = if tt1.abs() > tt2.abs() {
78        (tt1, tt2 + dtrd)
79    } else {
80        (tt1 + dtrd, tt2)
81    };
82    Ok(((tdb1, tdb2), 0))
83}
84
85// TT → UT1 (dt = TT−UT1 seconds).
86pub fn eraTtut1_safe(tt1: f64, tt2: f64, dt: f64) -> ErfaResult<((f64, f64), i32)> {
87    let dtd = dt / ERFA_DAYSEC;
88    let (ut11, ut12) = if tt1.abs() > tt2.abs() {
89        (tt1, tt2 - dtd)
90    } else {
91        (tt1 - dtd, tt2)
92    };
93    Ok(((ut11, ut12), 0))
94}
95
96// UT1 → TAI (dta = UT1−TAI seconds).
97pub fn eraUt1tai_safe(ut11: f64, ut12: f64, dta: f64) -> ErfaResult<((f64, f64), i32)> {
98    let dtad = dta / ERFA_DAYSEC;
99    let (tai1, tai2) = if ut11.abs() > ut12.abs() {
100        (ut11, ut12 - dtad)
101    } else {
102        (ut11 - dtad, ut12)
103    };
104    Ok(((tai1, tai2), 0))
105}
106
107// UT1 → TT (dt = TT−UT1 seconds).
108pub fn eraUt1tt_safe(ut11: f64, ut12: f64, dt: f64) -> ErfaResult<((f64, f64), i32)> {
109    let dtd = dt / ERFA_DAYSEC;
110    let (tt1, tt2) = if ut11.abs() > ut12.abs() {
111        (ut11, ut12 + dtd)
112    } else {
113        (ut11 + dtd, ut12)
114    };
115    Ok(((tt1, tt2), 0))
116}
117
118// eraUt1utc_safe  UT1 → UTC (handles leap-second wrinkles)
119// Returns: ((utc1, utc2), js) where js = 0 OK, +1 dubious year, -1 error
120pub fn eraUt1utc_safe(ut11: f64, ut12: f64, dut1: f64) -> ErfaResult<((f64, f64), i32)> {
121    // Arrange inputs big-first
122    let big1 = ut11.abs() >= ut12.abs();
123    let (u1, mut u2) = if big1 { (ut11, ut12) } else { (ut12, ut11) };
124
125    // Start with caller-supplied UT1-UTC.
126    let mut duts = dut1;
127
128    // Variables used inside the leap-second detection loop.
129    let d1 = u1; // first JD part (big)
130    let mut dats1: f64 = 0.0; // TAI-UTC at i = 1
131    let mut js: i32 = 0; // status from eraDat_safe
132
133    // Scan for possible leap-second day: from -1 to +3 days of UT1.
134    for i in -1..=3 {
135        let d2 = u2 + i as f64;
136
137        // Convert candidate day to calendar date.
138        // NOTE: fractional day '_fd' is intentionally unused in this path.
139        let ((iy, im, id), _fd, jcal) = eraJd2cal_safe(d1, d2)?;
140        if jcal != 0 {
141            return Ok(((0.0, 0.0), -1));
142        }
143
144        // TAI-UTC at 0h of this candidate day.
145        let (dats2, jdat) = eraDat_safe(iy, im, id, 0.0)?;
146        if jdat < 0 {
147            return Ok(((0.0, 0.0), -1));
148        }
149        js = jdat; // track dubious year flag if any
150
151        // Record first value (i = 1).
152        if i == -1 {
153            dats1 = dats2;
154        }
155
156        // Difference between consecutive days TAI-UTC.
157        let ddats = dats2 - dats1;
158
159        // Jump of ≥0.5 s implies leap-second boundary.
160        if ddats.abs() >= 0.5 {
161            // Ensure duts is the before value.
162            if ddats * duts >= 0.0 {
163                duts -= ddats;
164            }
165
166            // JD(UTC) of 0h UTC that ends with the leap second.
167            let ((us1, mut us2), jcz) = eraCal2jd_safe(iy, im, id)?;
168            if jcz != 0 {
169                return Ok(((0.0, 0.0), -1));
170            }
171            // Subtract 1 day then add current duts (now before value).
172            us2 = us2 - 1.0 + duts / ERFA_DAYSEC;
173
174            // How far into the current UTC day is the given UT1?
175            let mut du = u1 - us1;
176            du += u2 - us2;
177
178            if du > 0.0 {
179                // Fraction of the UTC day that has elapsed.
180                let fd = du * ERFA_DAYSEC / (ERFA_DAYSEC + ddats);
181                // Ramp UT1-UTC linearly during the leap second.
182                duts += ddats * if fd <= 1.0 { fd } else { 1.0 };
183            }
184
185            // Leap-second processing finished.
186            break;
187        }
188
189        // Prepare for next iteration.
190        dats1 = dats2;
191    }
192
193    // Subtract (possibly adjusted) UT1-UTC to obtain UTC.
194    u2 -= duts / ERFA_DAYSEC;
195
196    // Restore original part ordering for result.
197    let (utc1, utc2) = if big1 { (u1, u2) } else { (u2, u1) };
198    Ok(((utc1, utc2), js))
199}
200
201// UTC → TAI (with leap-second handling).
202pub fn eraUtctai_safe(utc1: f64, utc2: f64) -> ErfaResult<((f64, f64), i32)> {
203    let big1 = utc1.abs() >= utc2.abs();
204    let (u1, u2) = if big1 { (utc1, utc2) } else { (utc2, utc1) };
205
206    // Calendar for UTC
207    let ((iy, im, id), mut fd, jcal) = eraJd2cal_safe(u1, u2)?;
208    if jcal != 0 {
209        return Ok(((0.0, 0.0), jcal));
210    }
211
212    // TAI-UTC at 0h
213    let (dat0, mut j) = eraDat_safe(iy, im, id, 0.0)?;
214    if j < 0 {
215        return Ok(((0.0, 0.0), j));
216    }
217
218    let (dat12, j12) = eraDat_safe(iy, im, id, 0.5)?;
219    if j12 < 0 {
220        return Ok(((0.0, 0.0), j12));
221    }
222
223    // TAI-UTC at 24h (next day)
224    let ((iyt, imt, idt), _w, jcal2) = eraJd2cal_safe(u1 + 1.5, u2 - fd)?;
225    if jcal2 != 0 {
226        return Ok(((0.0, 0.0), jcal2));
227    }
228    let (dat24, j24) = eraDat_safe(iyt, imt, idt, 0.0)?;
229    if j24 < 0 {
230        return Ok(((0.0, 0.0), j24));
231    }
232    j = j24;
233
234    // Interpolate for any day duration changes
235    let dlod = 2.0 * (dat12 - dat0);
236    let dleap = dat24 - (dat0 + dlod);
237
238    fd *= (ERFA_DAYSEC + dleap) / ERFA_DAYSEC;
239    fd *= (ERFA_DAYSEC + dlod) / ERFA_DAYSEC;
240
241    // Build TAI parts
242    let ((z1, z2), jcz) = eraCal2jd_safe(iy, im, id)?;
243    if jcz != 0 {
244        return Ok(((0.0, 0.0), -1));
245    }
246    let a2 = z1 - u1 + z2 + fd + dat0 / ERFA_DAYSEC;
247
248    let (tai1, tai2) = if big1 { (u1, a2) } else { (a2, u1) };
249    Ok(((tai1, tai2), j))
250}
251
252// eraUtcut1_safe  UTC → UT1
253// Returns: ((ut11, ut12), js) where js=0 OK, +1 dubious year, -1 error
254pub fn eraUtcut1_safe(utc1: f64, utc2: f64, dut1: f64) -> ErfaResult<((f64, f64), i32)> {
255    // Date
256    let ((iy, im, id), _w, jcal) = eraJd2cal_safe(utc1, utc2)?;
257    if jcal != 0 {
258        return Ok(((0.0, 0.0), -1));
259    }
260
261    // TAI-UTC at 0h for the day
262    let (dat, mut js) = eraDat_safe(iy, im, id, 0.0)?;
263    if js < 0 {
264        return Ok(((0.0, 0.0), -1));
265    }
266
267    // UT1-TAI = (UT1-UTC) - (TAI-UTC)
268    let dta = dut1 - dat;
269
270    // UTC → TAI
271    let ((tai1, tai2), jw) = eraUtctai_safe(utc1, utc2)?;
272    if jw < 0 {
273        return Ok(((0.0, 0.0), -1));
274    } else if jw > 0 {
275        js = jw;
276    }
277
278    // TAI → UT1
279    let ((ut11, ut12), jt) = eraTaiut1_safe(tai1, tai2, dta)?;
280    if jt != 0 {
281        return Ok(((0.0, 0.0), -1));
282    }
283
284    Ok(((ut11, ut12), js))
285}