erfa_rust/
G19_safe.rs

1// G19
2//   icrs2g.c  → eraIcrs2g_safe
3//   ir.c      → eraIr_safe
4//   jd2cal.c  → eraJd2cal_safe
5//   jdcalf.c  → eraJdcalf_safe
6
7use crate::G1_safe::{eraAnp_safe, eraAnpm_safe};
8use crate::G28_safe::eraRxp_safe;
9use crate::G29_safe::eraS2c_safe;
10use crate::G7_safe::eraC2s_safe;
11use crate::H1_safe::ERFA_DNINT;
12
13pub type ErfaResult<T> = Result<T, ()>;
14
15//----------------------------------------------------------------------
16// G19/icrs2g.c → eraIcrs2g_safe
17//----------------------------------------------------------------------
18// Convert ICRS (dr, dd) to Galactic (dl, db), radians.
19pub fn eraIcrs2g_safe(dr: f64, dd: f64) -> ErfaResult<(f64, f64)> {
20    // ICRS→Galactic rotation matrix (row-major).
21    const R: [[f64; 3]; 3] = [
22        [
23            -0.054_875_560_416_215_368_492_398_900_454,
24            -0.873_437_090_234_885_048_760_383_168_409,
25            -0.483_835_015_548_713_226_831_774_175_116,
26        ],
27        [
28            0.494_109_427_875_583_673_525_222_371_358,
29            -0.444_829_629_960_011_178_146_614_061_616,
30            0.746_982_244_497_218_890_527_388_004_556,
31        ],
32        [
33            -0.867_666_149_019_004_701_181_616_534_570,
34            -0.198_076_373_431_201_528_180_486_091_412,
35            0.455_983_776_175_066_922_272_100_478_348,
36        ],
37    ];
38
39    let v1 = eraS2c_safe(dr, dd)?;
40    let v2 = eraRxp_safe(&R, &v1)?;
41    let (mut dl, mut db) = eraC2s_safe(&v2)?;
42    dl = eraAnp_safe(dl)?;
43    db = eraAnpm_safe(db)?;
44    Ok((dl, db))
45}
46
47//----------------------------------------------------------------------
48// G19/ir.c → eraIr_safe
49//----------------------------------------------------------------------
50// Set 3×3 matrix to identity.
51pub fn eraIr_safe(r: &mut [[f64; 3]; 3]) -> ErfaResult<()> {
52    for i in 0..3 {
53        for j in 0..3 {
54            r[i][j] = if i == j { 1.0 } else { 0.0 };
55        }
56    }
57    Ok(())
58}
59
60//----------------------------------------------------------------------
61// G19/jd2cal.c → eraJd2cal_safe
62//----------------------------------------------------------------------
63// Two-part JD to Gregorian (iy, im, id) and fractional day fd.
64// Returns ((iy, im, id), fd, j) with j=0 OK, -1 out of range.
65pub fn eraJd2cal_safe(dj1: f64, dj2: f64) -> ErfaResult<((i32, i32, i32), f64, i32)> {
66    const DJMIN: f64 = -68_569.5;
67    const DJMAX: f64 = 1e9;
68
69    // Range check.
70    let dj = dj1 + dj2;
71    if dj < DJMIN || dj > DJMAX {
72        return Ok(((0, 0, 0), 0.0, -1));
73    }
74
75    // Separate integer and fractional parts (compensated summation).
76    let mut jd = ERFA_DNINT(dj1) as i64;
77    let f1 = dj1 - jd as f64;
78    let d = ERFA_DNINT(dj2);
79    let f2 = dj2 - d;
80    jd += d as i64;
81
82    // Add 0.5 with compensation (KahanNeumaier).
83    let mut s = 0.5;
84    let mut cs = 0.0;
85    for x in [f1, f2] {
86        let t = s + x;
87        cs += if s.abs() >= x.abs() {
88            (s - t) + x
89        } else {
90            (x - t) + s
91        };
92        s = t;
93        if s >= 1.0 {
94            jd += 1;
95            s -= 1.0;
96        }
97    }
98    let mut f = s + cs;
99
100    // Handle negative fraction.
101    if f < 0.0 {
102        f += 1.0;
103        jd -= 1;
104    }
105
106    // Handle rounding up to 1.0.
107    if (f - 1.0) >= -f64::EPSILON / 4.0 {
108        f -= 1.0;
109        jd += 1;
110        if f < 0.0 {
111            f = 0.0;
112        }
113    }
114
115    // Gregorian calendar from integer JD (Fliegel & Van Flandern).
116    let mut l = jd + 68_569;
117    let n = (4 * l) / 146_097;
118    l -= (146_097 * n + 3) / 4;
119    let i = (4_000 * (l + 1)) / 1_461_001;
120    l -= (1_461 * i) / 4 - 31;
121    let k = (80 * l) / 2_447;
122
123    let id = (l - (2_447 * k) / 80) as i32;
124    l = k / 11;
125    let im = (k + 2 - 12 * l) as i32;
126    let iy = (100 * (n - 49) + i + l) as i32;
127
128    Ok(((iy, im, id), f, 0))
129}
130
131//----------------------------------------------------------------------
132// G19/jdcalf.c → eraJdcalf_safe
133//----------------------------------------------------------------------
134// JD to calendar with ndp decimals; returns [iy, im, id, frac×10^ndp] and status.
135pub fn eraJdcalf_safe(ndp: i32, dj1: f64, dj2: f64) -> ErfaResult<([i32; 4], i32)> {
136    let (mut jstat, denom) = if (0..=9).contains(&ndp) {
137        (0, 10_f64.powi(ndp))
138    } else {
139        (1, 1.0)
140    };
141
142    // Arrange parts so |d1| ≥ |d2|.
143    let (mut d1, d2) = if dj1.abs() >= dj2.abs() {
144        (dj1, dj2)
145    } else {
146        (dj2, dj1)
147    };
148
149    // Realign to midnight.
150    d1 -= 0.5;
151
152    // Separate day & fraction.
153    let mut d = ERFA_DNINT(d1);
154    let f1 = d1 - d;
155    let mut djd = d;
156    d = ERFA_DNINT(d2);
157    let f2 = d2 - d;
158    djd += d;
159    d = ERFA_DNINT(f1 + f2);
160    let mut f = (f1 - d) + f2;
161    if f < 0.0 {
162        f += 1.0;
163        d -= 1.0;
164    }
165    djd += d;
166
167    // Round fraction to requested decimals.
168    let rf = ERFA_DNINT(f * denom) / denom;
169
170    // Re-align to noon.
171    djd += 0.5;
172
173    // Convert to calendar date.
174    let ((year, month, day), frac, js) = eraJd2cal_safe(djd, rf)?;
175    let mut out = [0i32; 4];
176    if js == 0 {
177        let frac_i = ERFA_DNINT(frac * denom) as i32;
178        out[0] = year;
179        out[1] = month;
180        out[2] = day;
181        out[3] = frac_i;
182    } else {
183        jstat = js;
184    }
185
186    Ok((out, jstat))
187}
erfa_rust/G19_safe.rs

erfa_rust/
G19_safe.rs
