erfa_rust/
G34_safe.rs

1// G34
2//   xy06.c  → eraXy06_safe
3
4use crate::G15_safe::{
5    eraFad03_safe, eraFae03_safe, eraFaf03_safe, eraFaju03_safe, eraFal03_safe, eraFalp03_safe,
6    eraFama03_safe, eraFame03_safe, eraFane03_safe, eraFaom03_safe, eraFapa03_safe, eraFasa03_safe,
7    eraFaur03_safe, eraFave03_safe,
8};
9use crate::H1_safe::{ERFA_DAS2R, ERFA_DJ00, ERFA_DJC};
10
11pub type ErfaResult<T> = Result<T, ()>;
12
13#[path = "data/G34_safe/MFALS.rs"]
14mod mfals_mod;
15use mfals_mod::MFALS;
16#[path = "data/G34_safe/MFAPL.rs"]
17mod mfapl_mod;
18use mfapl_mod::MFAPL;
19#[path = "data/G34_safe/NC.rs"]
20mod nc_mod;
21use nc_mod::NC;
22#[path = "data/G34_safe/AMPL.rs"]
23mod ampl_mod;
24use ampl_mod::AMPL as A;
25
26// Maximum power of T in the X,Y precession-nutation polynomials
27const MAXPT: usize = 5;
28
29// Polynomial coefficients (arcsec).  First row = X, second = Y
30static XYP: [[f64; MAXPT + 1]; 2] = [
31    [
32        -0.016_617,
33        2004.191_898,
34        -0.429_782_9,
35        -0.198_618_34,
36        0.000_007_578,
37        0.000_005_928_5,
38    ],
39    [
40        -0.006_951,
41        -0.025_896,
42        -22.407_274_7,
43        0.001_900_59,
44        0.001_112_526,
45        0.000_000_135_8,
46    ],
47];
48
49/* Amplitude usage: X/Y, sin/cos, power-of-T -- length=20 each */
50static JAXY: [usize; 20] = [0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1];
51static JASC: [usize; 20] = [0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0];
52static JAPT: [usize; 20] = [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4];
53
54const NFLS: usize = MFALS.len();
55const NFPL: usize = MFAPL.len();
56const NA: usize = A.len();
57
58// CIP X,Y using IAU-2006/2000A series (truncated data).
59pub fn eraXy06_safe(date1: f64, date2: f64) -> ErfaResult<(f64, f64)> {
60    // Time interval in Julian-centuries since J2000.0
61    let t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
62
63    // Powers of T
64    let mut pt = [0.0_f64; MAXPT + 1];
65    pt[0] = 1.0;
66    for j in 1..=MAXPT {
67        pt[j] = pt[j - 1] * t;
68    }
69
70    // Accumulators: polynomial (pr), luni-solar (ls), planetary (pl)
71    let mut xypr = [0.0_f64; 2];
72    let mut xyls = [0.0_f64; 2];
73    let mut xypl = [0.0_f64; 2];
74
75    // Fundamental arguments (IERS 2003)
76    let mut fa = [0.0_f64; 14];
77    fa[0] = eraFal03_safe(t)?; // l
78    fa[1] = eraFalp03_safe(t)?; // l'
79    fa[2] = eraFaf03_safe(t)?; // F
80    fa[3] = eraFad03_safe(t)?; // D
81    fa[4] = eraFaom03_safe(t)?; // Om
82    fa[5] = eraFame03_safe(t)?; // Mercury
83    fa[6] = eraFave03_safe(t)?; // Venus
84    fa[7] = eraFae03_safe(t)?; // Earth
85    fa[8] = eraFama03_safe(t)?; // Mars
86    fa[9] = eraFaju03_safe(t)?; // Jupiter
87    fa[10] = eraFasa03_safe(t)?; // Saturn
88    fa[11] = eraFaur03_safe(t)?; // Uranus
89    fa[12] = eraFane03_safe(t)?; // Neptune
90    fa[13] = eraFapa03_safe(t)?; // precession
91
92    // Polynomial part (reverse iteration)
93    for jxy in 0..2 {
94        for j in (0..=MAXPT).rev() {
95            xypr[jxy] += XYP[jxy][j] * pt[j];
96        }
97    }
98
99    // Helper: accumulate one term into provided target array
100    let accumulate = |target: &mut [f64; 2],
101                      jxy: usize,
102                      jsc: usize,
103                      jpt: usize,
104                      amp_index: usize,
105                      sc: [f64; 2]| {
106        if amp_index < NA {
107            target[jxy] += A[amp_index] * sc[jsc] * pt[jpt];
108        }
109    };
110
111    // Planetary nutation terms
112    let mut ialast: isize = NA as isize; // 1-based like C code
113    for ifreq in (0..NFPL).rev() {
114        // Argument = Σ m_i * fa_i
115        let mut arg = 0.0_f64;
116        for i in 0..14 {
117            let m = MFAPL[ifreq][i];
118            if m != 0 {
119                arg += (m as f64) * fa[i];
120            }
121        }
122        let sc = [arg.sin(), arg.cos()];
123
124        // Iterate over amplitudes for this frequency
125        let ia = NC[ifreq + NFLS] as isize;
126        for i in (ia..=ialast).rev() {
127            let j = (i - ia) as usize;
128            if j < 20 {
129                accumulate(&mut xypl, JAXY[j], JASC[j], JAPT[j], (i - 1) as usize, sc);
130            }
131        }
132        ialast = ia - 1;
133    }
134
135    // Luni-solar nutation terms
136    for ifreq in (0..NFLS).rev() {
137        // Argument = Σ m_i * fa_i
138        let mut arg = 0.0_f64;
139        for i in 0..5 {
140            let m = MFALS[ifreq][i];
141            if m != 0 {
142                arg += (m as f64) * fa[i];
143            }
144        }
145        let sc = [arg.sin(), arg.cos()];
146
147        // Iterate over amplitudes for this frequency
148        let ia = NC[ifreq] as isize;
149        for i in (ia..=ialast).rev() {
150            let j = (i - ia) as usize;
151            if j < 20 {
152                accumulate(&mut xyls, JAXY[j], JASC[j], JAPT[j], (i - 1) as usize, sc);
153            }
154        }
155        ialast = ia - 1;
156    }
157
158    // Result: CIP unit-vector components (radians)
159    let x = ERFA_DAS2R * (xypr[0] + (xyls[0] + xypl[0]) / 1.0e6);
160    let y = ERFA_DAS2R * (xypr[1] + (xyls[1] + xypl[1]) / 1.0e6);
161
162    Ok((x, y))
163}
erfa_rust/G34_safe.rs

erfa_rust/
G34_safe.rs
