erfa_rust/
G10_safe.rs

1// G10
2//   dtdb.c  → eraDtdb_safe
3//   dtf2d.c → eraDtf2d_safe
4
5use crate::G19_safe::eraJd2cal_safe;
6use crate::G8_safe::eraCal2jd_safe;
7use crate::G9_safe::eraDat_safe;
8use crate::H1_safe::{ERFA_D2PI, ERFA_DAYSEC, ERFA_DD2R, ERFA_DJ00, ERFA_DJM};
9
10#[path = "data/G10_safe/FAIRHD.rs"]
11mod fairhd_mod;
12use fairhd_mod::FAIRHD;
13
14pub type ErfaResult<T> = Result<T, ()>;
15
16#[inline]
17fn fmod(a: f64, b: f64) -> f64 {
18    // Emulate C fmod with correct sign for negatives
19    a - (a / b).trunc() * b
20}
21
22// Geocentric model for TDB − TT (seconds).
23pub fn eraDtdb_safe(
24    date1: f64,
25    date2: f64,
26    ut: f64,
27    elong: f64,
28    u: f64,
29    v: f64,
30) -> ErfaResult<f64> {
31    // Time since J2000.0 in Julian millennia
32    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJM;
33
34    // Topocentric part (Moyer/Murray)
35    let tsol = fmod(ut, 1.0) * ERFA_D2PI + elong;
36
37    // Fundamental arguments  Simon et al. 1994
38    let w = t / 3600.0;
39
40    // Mean longitudes and anomalies (radians)
41    let elsun = fmod(280.466_456_83 + 1_296_027_711.034_29 * w, 360.0) * ERFA_DD2R;
42    let emsun = fmod(357.529_109_18 + 1_295_965_810.481_00 * w, 360.0) * ERFA_DD2R;
43    let d = fmod(297.850_195_47 + 16_029_616_012.090_00 * w, 360.0) * ERFA_DD2R;
44    let elj = fmod(34.351_518_74 + 109_306_899.894_53 * w, 360.0) * ERFA_DD2R;
45    let els = fmod(50.077_444_30 + 44_046_398.470_38 * w, 360.0) * ERFA_DD2R;
46
47    // Moyer + Murray topocentric terms
48    let wt = 0.000_29e-10 * u * (tsol + elsun - els).sin()
49        + 0.001_00e-10 * u * (tsol - 2.0 * emsun).sin()
50        + 0.001_33e-10 * u * (tsol - d).sin()
51        + 0.001_33e-10 * u * (tsol + elsun - elj).sin()
52        - 0.002_29e-10 * u * (tsol + 2.0 * elsun + emsun).sin()
53        - 0.022_00e-10 * v * (elsun + emsun).cos()
54        + 0.053_12e-10 * u * (tsol - emsun).sin()
55        - 0.136_77e-10 * u * (tsol + 2.0 * elsun).sin()
56        - 1.318_40e-10 * v * elsun.cos()
57        + 3.176_79e-10 * u * tsol.sin();
58
59    // Fairhead & Bretagnon (geocentric): T^0 .. T^4 harmonic sums
60    let mut w0 = 0.0;
61    for j in (0..=473).rev() {
62        let row = FAIRHD[j];
63        w0 += row[0] * (row[1] * t + row[2]).sin();
64    }
65    let mut w1 = 0.0;
66    for j in (474..=678).rev() {
67        let row = FAIRHD[j];
68        w1 += row[0] * (row[1] * t + row[2]).sin();
69    }
70    let mut w2 = 0.0;
71    for j in (679..=763).rev() {
72        let row = FAIRHD[j];
73        w2 += row[0] * (row[1] * t + row[2]).sin();
74    }
75    let mut w3 = 0.0;
76    for j in (764..=783).rev() {
77        let row = FAIRHD[j];
78        w3 += row[0] * (row[1] * t + row[2]).sin();
79    }
80    let mut w4 = 0.0;
81    for j in (784..=786).rev() {
82        let row = FAIRHD[j];
83        w4 += row[0] * (row[1] * t + row[2]).sin();
84    }
85
86    // Combine powers of T
87    let wf = t * (t * (t * (t * w4 + w3) + w2) + w1) + w0;
88
89    // JPL-mass adjustments
90    let wj = 0.000_65e-6 * (6_069.776_754 * t + 4.021_194).sin()
91        + 0.000_33e-6 * (213.299_095 * t + 5.543_132).sin()
92        - 0.001_96e-6 * (6_208.294_251 * t + 5.696_701).sin()
93        - 0.001_73e-6 * (74.781_599 * t + 2.435_900).sin()
94        + 0.036_38e-6 * t * t;
95
96    // Final result: TDB − TT (seconds)
97    Ok(wt + wf + wj)
98}
99
100// Calendar + clock to two-part Julian Date/quasi-JD.
101pub fn eraDtf2d_safe(
102    scale: &str,
103    iy: i32,
104    im: i32,
105    id: i32,
106    ihr: i32,
107    imn: i32,
108    sec: f64,
109) -> Result<((f64, f64), i32), i32> {
110    // Convert calendar date to JD (0h today)
111    let ((djm0, djm), j_cal) = eraCal2jd_safe(iy, im, id).map_err(|_| -9999)?;
112    if j_cal != 0 {
113        return Err(j_cal);
114    }
115    let dj = djm0 + djm;
116
117    // Default day length and final-minute seconds
118    let mut day = ERFA_DAYSEC;
119    let mut seclim = 60.0_f64;
120
121    // Leap-second handling for UTC only (exact match)
122    if scale == "UTC" {
123        // TAI−UTC at 0h today
124        let (dat0, j0) = eraDat_safe(iy, im, id, 0.0).map_err(|_| -9999)?;
125        if j0 < 0 {
126            return Err(j0);
127        }
128
129        // TAI−UTC at 12h today
130        let (dat12, j12) = eraDat_safe(iy, im, id, 0.5).map_err(|_| -9999)?;
131        if j12 < 0 {
132            return Err(j12);
133        }
134
135        // TAI−UTC at 0h tomorrow
136        let ((iy2, im2, id2), _w, j_next) = eraJd2cal_safe(dj + 1.5, 0.0).map_err(|_| -9999)?;
137        if j_next != 0 {
138            return Err(j_next);
139        }
140        let (dat24, j24) = eraDat_safe(iy2, im2, id2, 0.0).map_err(|_| -9999)?;
141        if j24 < 0 {
142            return Err(j24);
143        }
144
145        // Leap-second increment for today
146        let dleap = dat24 - (2.0 * dat12 - dat0);
147
148        // Adjust day and final-minute seconds on leap-second days
149        day += dleap;
150        if ihr == 23 && imn == 59 {
151            seclim += dleap;
152        }
153    }
154
155    // Validate clock fields and produce warning if time after end-of-day
156    let mut js_warn = 0;
157    if (0..=23).contains(&ihr) {
158        if (0..=59).contains(&imn) {
159            if sec >= 0.0 {
160                if sec >= seclim {
161                    js_warn += 2;
162                }
163            } else {
164                return Err(-6);
165            }
166        } else {
167            return Err(-5);
168        }
169    } else {
170        return Err(-4);
171    }
172
173    // Fraction of the day
174    let time = (60.0 * (60.0 * ihr as f64 + imn as f64) + sec) / day;
175
176    // Return ((d1, d2), j) with j = 0 or +2
177    Ok(((dj, time), js_warn))
178}
erfa_rust/G10_safe.rs

erfa_rust/
G10_safe.rs
