stem_branch/
lib.rs

1//! Native Rust port of stem-branch's solar ephemeris core.
2//!
3//! Computes the Sun's geocentric ecliptic state from the full VSOP87D Earth
4//! series plus a JPL DE441-fitted correction polynomial and IAU2000B nutation.
5//! Input is a Julian Ephemeris Day in Terrestrial Time (JDE/TT); the UTC<->TT /
6//! ΔT boundary is the caller's responsibility.
7//!
8//! This is a faithful translation of the TypeScript `solarEclipticState`; the
9//! VSOP87D and IAU2000B coefficient tables are generated from the same source
10//! by `scripts/gen-rust-solar-data.mjs`.
11
12#![forbid(unsafe_code)]
13
14mod delta_t;
15mod elpmpp02_data;
16mod julian;
17mod lunisolar;
18mod moon;
19mod new_moon;
20mod nutation_data;
21mod solar_terms;
22mod vsop87d_earth;
23
24pub use delta_t::delta_t_for_year;
25pub use julian::{jd_from_ymd, ymd_from_jd};
26pub use lunisolar::{
27    gregorian_to_lunisolar, lunar_months_for_year, lunar_new_year, CivilDate, LunarMonth,
28    LunisolarDate,
29};
30pub use moon::{moon_position, MoonState};
31pub use new_moon::{find_new_moons_in_range, new_moon_jde};
32pub use solar_terms::{find_solar_term_moment, SOLAR_TERM_LONGITUDES};
33
34use core::f64::consts::{PI, TAU};
35use nutation_data::{NUT_COEFFS, NUT_OBLIQ};
36use vsop87d_earth::{EARTH_L, EARTH_R};
37
38/// Radians per arcsecond (π / 180 / 3600).
39pub(crate) const ARCSEC_TO_RAD: f64 = PI / 180.0 / 3600.0;
40/// Degrees per radian.
41pub(crate) const RAD_TO_DEG: f64 = 180.0 / PI;
42
43/// Geocentric solar ecliptic state at a Julian Ephemeris Day (TT).
44#[derive(Debug, Clone, Copy, PartialEq)]
45pub struct SolarState {
46    /// Geocentric ecliptic longitude, mean equinox of date (no nutation/
47    /// aberration), in degrees `[0, 360)`.
48    pub true_longitude_degrees: f64,
49    /// Apparent geocentric ecliptic longitude — `true` plus IAU2000B nutation
50    /// in longitude and aberration (true equinox of date), in degrees `[0, 360)`.
51    pub apparent_longitude_degrees: f64,
52    /// Sun–Earth distance in astronomical units.
53    pub radius_au: f64,
54}
55
56/// Compute the Sun's geocentric ecliptic state at the given Julian Ephemeris
57/// Day in Terrestrial Time.
58pub fn solar_ecliptic_state(jde_tt: f64) -> SolarState {
59    let tau = (jde_tt - 2451545.0) / 365250.0; // Julian millennia from J2000 (TT)
60    let t = (jde_tt - 2451545.0) / 36525.0; // Julian centuries from J2000 (TT)
61
62    // VSOP87D heliocentric longitude (rad) and radius (AU).
63    let mut lon = eval_vsop_series(EARTH_L, tau);
64    let r = eval_vsop_series(EARTH_R, tau);
65
66    // DE441-fitted even-polynomial correction (arcseconds -> radians). The
67    // literal 206264.806 mirrors the TypeScript source.
68    let tau2 = tau * tau;
69    lon += (-0.106674 - 0.616597 * tau2 + 0.315446 * tau2 * tau2 - 0.050315 * tau2 * tau2 * tau2)
70        / 206264.806;
71
72    // Heliocentric -> geocentric (+180°): mean equinox of date before nutation
73    // and aberration.
74    let geo_true = lon + PI;
75
76    // Apparent place: + IAU2000B nutation in longitude + aberration.
77    let (dl, dlp, df, dd, dom) = delaunay_args(t);
78    let dpsi = nutation_dpsi(dl, dlp, df, dd, dom, t);
79    let apparent = geo_true + dpsi * ARCSEC_TO_RAD + (-20.4898 / r) * ARCSEC_TO_RAD;
80
81    SolarState {
82        true_longitude_degrees: normalize_radians(geo_true) * RAD_TO_DEG,
83        apparent_longitude_degrees: normalize_radians(apparent) * RAD_TO_DEG,
84        radius_au: r,
85    }
86}
87
88/// Evaluate a VSOP87 series: a sum over powers of `tau`, each power weighting a
89/// sum of `A * cos(B + C * tau)` terms.
90fn eval_vsop_series(series: &[&[[f64; 3]]], tau: f64) -> f64 {
91    let mut result = 0.0;
92    let mut tau_pow = 1.0;
93    for terms in series {
94        let mut sum = 0.0;
95        for term in *terms {
96            sum += term[0] * (term[1] + term[2] * tau).cos();
97        }
98        result += sum * tau_pow;
99        tau_pow *= tau;
100    }
101    result
102}
103
104/// The five Delaunay fundamental arguments (radians) at Julian century `t` (TT).
105/// Returns `(l, l', F, D, Ω)`. Source: IERS Conventions (2010), Table 5.2a.
106pub(crate) fn delaunay_args(t: f64) -> (f64, f64, f64, f64, f64) {
107    let t2 = t * t;
108    let t3 = t2 * t;
109    let t4 = t3 * t;
110    // `% 1296000.0` (arcseconds in 360°) reduces each argument; the f64 `%`
111    // operator matches JavaScript's truncated remainder for a faithful port.
112    let l = ((485868.249036 + 1717915923.2178 * t + 31.8792 * t2 + 0.051635 * t3
113        - 0.00024470 * t4)
114        % 1296000.0)
115        * ARCSEC_TO_RAD;
116    let lp = ((1287104.79305 + 129596581.0481 * t - 0.5532 * t2 + 0.000136 * t3 - 0.00001149 * t4)
117        % 1296000.0)
118        * ARCSEC_TO_RAD;
119    let f = ((335779.526232 + 1739527262.8478 * t - 12.7512 * t2 - 0.001037 * t3
120        + 0.00000417 * t4)
121        % 1296000.0)
122        * ARCSEC_TO_RAD;
123    let d = ((1072260.70369 + 1602961601.2090 * t - 6.3706 * t2 + 0.006593 * t3 - 0.00003169 * t4)
124        % 1296000.0)
125        * ARCSEC_TO_RAD;
126    let om = ((450160.398036 - 6962890.5431 * t + 7.4722 * t2 + 0.007702 * t3 - 0.00005939 * t4)
127        % 1296000.0)
128        * ARCSEC_TO_RAD;
129    (l, lp, f, d, om)
130}
131
132/// Nutation in longitude (Δψ) in arcseconds, IAU2000B (77 lunisolar terms).
133pub(crate) fn nutation_dpsi(l: f64, lp: f64, f: f64, d: f64, om: f64, t: f64) -> f64 {
134    let mut dpsi = 0.0;
135    for row in NUT_COEFFS {
136        let arg = row[0] * l + row[1] * lp + row[2] * f + row[3] * d + row[4] * om;
137        dpsi += (row[5] + row[6] * t) * arg.sin();
138    }
139    // 0.1 microarcseconds -> arcseconds.
140    dpsi / 1e7
141}
142
143/// Nutation in obliquity (Δε) in arcseconds, IAU2000B (parallel to `NUT_COEFFS`).
144pub(crate) fn nutation_deps(l: f64, lp: f64, f: f64, d: f64, om: f64, t: f64) -> f64 {
145    let mut deps = 0.0;
146    for (row, obliq) in NUT_COEFFS.iter().zip(NUT_OBLIQ) {
147        let arg = row[0] * l + row[1] * lp + row[2] * f + row[3] * d + row[4] * om;
148        deps += (obliq[0] + obliq[1] * t) * arg.cos();
149    }
150    // 0.1 microarcseconds -> arcseconds.
151    deps / 1e7
152}
153
154/// Normalize an angle in radians to `[0, 2π)`.
155pub(crate) fn normalize_radians(rad: f64) -> f64 {
156    ((rad % TAU) + TAU) % TAU
157}
stem_branch/lib.rs

stem_branch/
lib.rs
