erfa_rust/
G23_safe.rs

1// G23
2//   nut00b.c → eraNut00b_safe
3//   nut06a.c → eraNut06a_safe
4//   nut80.c  → eraNut80_safe
5//   nutm80.c → eraNutm80_safe
6//   obl06.c  → eraObl06_safe
7//   obl80.c  → eraObl80_safe
8
9use crate::G1_safe::eraAnpm_safe;
10use crate::G21_safe::eraNumat_safe;
11use crate::G22_safe::eraNut00a_safe;
12use crate::H1_safe::{ERFA_D2PI, ERFA_DAS2R, ERFA_DJ00, ERFA_DJC, ERFA_DMAS2R, ERFA_TURNAS};
13
14#[path = "data/G23_safe/X00B.rs"]
15mod x00b_mod;
16use x00b_mod::X00B;
17
18#[path = "data/G23_safe/X80.rs"]
19mod x80_mod;
20use x80_mod::X80;
21
22pub type ErfaResult<T> = Result<T, ()>;
23
24// Helper: positive modulo for angles.
25#[inline]
26fn fmod_pos(a: f64, b: f64) -> f64 {
27    a.rem_euclid(b)
28}
29
30
31 
32//  eraNut00b_safe   IAU-2000B nutation
33
34/* 77 luni-solar */
35#[derive(Clone, Copy)]
36struct Term00B {
37    nl: i32,
38    nlp: i32,
39    nf: i32,
40    nd: i32,
41    nom: i32,
42    ps: f64,
43    pst: f64,
44    pc: f64,
45    ec: f64,
46    ect: f64,
47    es: f64,
48}
49
50const U2R_00B: f64 = ERFA_DAS2R / 1e7; /* 0.1 µas → rad */
51const DPPLAN: f64 = -0.135 * ERFA_DMAS2R;
52const DEPLAN: f64 = 0.388 * ERFA_DMAS2R;
53
54pub fn eraNut00b_safe(date1: f64, date2: f64) -> ErfaResult<(f64, f64)> {
55    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
56
57    // Fundamental arguments (arcsec -> rad for linear terms).
58    let el = fmod_pos(485_868.249_036 + 1_717_915_923.2178 * t, ERFA_TURNAS) * ERFA_DAS2R;
59    let elp = fmod_pos(1_287_104.793_05 + 129_596_581.0481 * t, ERFA_TURNAS) * ERFA_DAS2R;
60    let f = fmod_pos(335_779.526_232 + 1_739_527_262.8478 * t, ERFA_TURNAS) * ERFA_DAS2R;
61    let d = fmod_pos(1_072_260.703_69 + 1_602_961_601.2090 * t, ERFA_TURNAS) * ERFA_DAS2R;
62    let om = fmod_pos(450_160.398_036 - 6_962_890.5431 * t, ERFA_TURNAS) * ERFA_DAS2R;
63
64    // Accumulate series.
65    let mut dp = 0.0f64;
66    let mut de = 0.0f64;
67    for term in X00B.iter().rev() {
68        let arg = fmod_pos(
69            (term.nl as f64) * el
70                + (term.nlp as f64) * elp
71                + (term.nf as f64) * f
72                + (term.nd as f64) * d
73                + (term.nom as f64) * om,
74            ERFA_D2PI,
75        );
76        let sarg = arg.sin();
77        let carg = arg.cos();
78        dp += (term.ps + term.pst * t) * sarg + term.pc * carg;
79        de += (term.ec + term.ect * t) * carg + term.es * sarg;
80    }
81
82    let dpsi = dp * U2R_00B + DPPLAN;
83    let deps = de * U2R_00B + DEPLAN;
84    Ok((dpsi, deps))
85}
86
87
88//  eraNut06a_safe   2000A nutation + 2006
89#[derive(Clone, Copy)]
90struct Term80 {
91    nl: i32,
92    nlp: i32,
93    nf: i32,
94    nd: i32,
95    nom: i32,
96    sp: f64,
97    spt: f64,
98    ce: f64,
99    cet: f64,
100}
101
102pub fn eraNut06a_safe(date1: f64, date2: f64) -> ErfaResult<(f64, f64)> {
103    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
104    let fj2 = -2.7774e-6 * t;
105
106    let (dp0, de0) = eraNut00a_safe(date1, date2)?;
107    let dpsi = dp0 + dp0 * (0.4697e-6 + fj2);
108    let deps = de0 + de0 * fj2;
109    Ok((dpsi, deps))
110}
111
112// IAU 1980 nutation: returns (Δψ, Δε) in radians.
113pub fn eraNut80_safe(date1: f64, date2: f64) -> ErfaResult<(f64, f64)> {
114    const U2R_80: f64 = ERFA_DAS2R / 1e4; // 0.1 mas → rad
115    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
116
117    // Delaunay arguments, normalized to ±π using eraAnpm_safe.
118    let el = eraAnpm_safe(
119        (485_866.733 + (715_922.633 + (31.310 + 0.064 * t) * t) * t) * ERFA_DAS2R
120            + fmod_pos(1325.0 * t, 1.0) * ERFA_D2PI,
121    )?;
122    let elp = eraAnpm_safe(
123        (1_287_099.804 + (1_292_581.224 + (-0.577 - 0.012 * t) * t) * t) * ERFA_DAS2R
124            + fmod_pos(99.0 * t, 1.0) * ERFA_D2PI,
125    )?;
126    let f = eraAnpm_safe(
127        (335_778.877 + (295_263.137 + (-13.257 + 0.011 * t) * t) * t) * ERFA_DAS2R
128            + fmod_pos(1342.0 * t, 1.0) * ERFA_D2PI,
129    )?;
130    let d = eraAnpm_safe(
131        (1_072_261.307 + (1_105_601.328 + (-6.891 + 0.019 * t) * t) * t) * ERFA_DAS2R
132            + fmod_pos(1236.0 * t, 1.0) * ERFA_D2PI,
133    )?;
134    let om = eraAnpm_safe(
135        (450_160.280 + (-482_890.539 + (7.455 + 0.008 * t) * t) * t) * ERFA_DAS2R
136            + fmod_pos(-5.0 * t, 1.0) * ERFA_D2PI,
137    )?;
138
139    // Accumulate series.
140    let mut dp = 0.0f64;
141    let mut de = 0.0f64;
142    for term in X80.iter().rev() {
143        let arg = (term.nl as f64) * el
144            + (term.nlp as f64) * elp
145            + (term.nf as f64) * f
146            + (term.nd as f64) * d
147            + (term.nom as f64) * om;
148        let s = term.sp + term.spt * t;
149        let c = term.ce + term.cet * t;
150        if s != 0.0 {
151            dp += s * arg.sin();
152        }
153        if c != 0.0 {
154            de += c * arg.cos();
155        }
156    }
157
158    Ok((dp * U2R_80, de * U2R_80))
159}
160
161// Nutation matrix (1980 model): builds 3×3 matrix from Δψ, Δε and mean obliquity.
162pub fn eraNutm80_safe(date1: f64, date2: f64) -> ErfaResult<[[f64; 3]; 3]> {
163    let (dp, de) = eraNut80_safe(date1, date2)?;
164    let epsa = eraObl80_safe(date1, date2)?;
165    let rmatn = eraNumat_safe(epsa, dp, de)?;
166    Ok(rmatn)
167}
168
169// Mean obliquity of the ecliptic, IAU 2006, radians.
170pub fn eraObl06_safe(date1: f64, date2: f64) -> ErfaResult<f64> {
171    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
172    let eps = (84_381.406
173        + (-46.836_769
174            + (-0.000_1831 + (0.002_003_40 + (-0.000_000_576 + (-0.000_000_0434) * t) * t) * t)
175                * t)
176            * t)
177        * ERFA_DAS2R;
178    Ok(eps)
179}
180
181// Mean obliquity of the ecliptic, IAU 1980, radians.
182pub fn eraObl80_safe(date1: f64, date2: f64) -> ErfaResult<f64> {
183    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
184    let eps = ERFA_DAS2R * (84_381.448 + (-46.8150 + (-0.00059 + 0.001_813 * t) * t) * t);
185    Ok(eps)
186}
erfa_rust/G23_safe.rs

erfa_rust/
G23_safe.rs
