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 names;
20mod new_moon;
21mod nutation_data;
22mod solar_terms;
23mod vsop87d_earth;
24
25pub use delta_t::delta_t_for_year;
26pub use julian::{jd_from_ymd, ymd_from_jd};
27pub use lunisolar::{
28    gregorian_to_lunisolar, gregorian_to_lunisolar_with, lunar_months_for_year, lunar_new_year,
29    CivilDate, LunarMonth, LunisolarDate,
30};
31pub use moon::{moon_position, MoonState};
32pub use names::{EARTHLY_BRANCHES, HEAVENLY_STEMS};
33pub use new_moon::{find_new_moons_in_range, new_moon_jde};
34pub use solar_terms::{
35    find_solar_term_moment, solar_term_for_longitude, SOLAR_TERM_LONGITUDES, SOLAR_TERM_NAMES,
36};
37
38use core::f64::consts::{PI, TAU};
39use nutation_data::{NUT_COEFFS, NUT_OBLIQ};
40use vsop87d_earth::{EARTH_L, EARTH_R};
41
42/// Radians per arcsecond (π / 180 / 3600).
43pub(crate) const ARCSEC_TO_RAD: f64 = PI / 180.0 / 3600.0;
44/// Degrees per radian.
45pub(crate) const RAD_TO_DEG: f64 = 180.0 / PI;
46
47/// Geocentric solar ecliptic state at a Julian Ephemeris Day (TT).
48#[derive(Debug, Clone, Copy, PartialEq)]
49pub struct SolarState {
50    /// Geocentric ecliptic longitude, mean equinox of date (no nutation/
51    /// aberration), in degrees `[0, 360)`.
52    pub true_longitude_degrees: f64,
53    /// Apparent geocentric ecliptic longitude — `true` plus IAU2000B nutation
54    /// in longitude and aberration (true equinox of date), in degrees `[0, 360)`.
55    pub apparent_longitude_degrees: f64,
56    /// Sun–Earth distance in astronomical units.
57    pub radius_au: f64,
58}
59
60/// Compute the Sun's geocentric ecliptic state at the given Julian Ephemeris
61/// Day in Terrestrial Time.
62pub fn solar_ecliptic_state(jde_tt: f64) -> SolarState {
63    let tau = (jde_tt - 2451545.0) / 365250.0; // Julian millennia from J2000 (TT)
64    let t = (jde_tt - 2451545.0) / 36525.0; // Julian centuries from J2000 (TT)
65
66    // VSOP87D heliocentric longitude (rad) and radius (AU).
67    let mut lon = eval_vsop_series(EARTH_L, tau);
68    let r = eval_vsop_series(EARTH_R, tau);
69
70    // DE441-fitted even-polynomial correction (arcseconds -> radians). The
71    // literal 206264.806 mirrors the TypeScript source.
72    let tau2 = tau * tau;
73    lon += (-0.106674 - 0.616597 * tau2 + 0.315446 * tau2 * tau2 - 0.050315 * tau2 * tau2 * tau2)
74        / 206264.806;
75
76    // Heliocentric -> geocentric (+180°): mean equinox of date before nutation
77    // and aberration.
78    let geo_true = lon + PI;
79
80    // Apparent place: + IAU2000B nutation in longitude + aberration.
81    let (dl, dlp, df, dd, dom) = delaunay_args(t);
82    let dpsi = nutation_dpsi(dl, dlp, df, dd, dom, t);
83    let apparent = geo_true + dpsi * ARCSEC_TO_RAD + (-20.4898 / r) * ARCSEC_TO_RAD;
84
85    SolarState {
86        true_longitude_degrees: normalize_radians(geo_true) * RAD_TO_DEG,
87        apparent_longitude_degrees: normalize_radians(apparent) * RAD_TO_DEG,
88        radius_au: r,
89    }
90}
91
92/// Astronomical unit in kilometres (IAU 2012 definition).
93pub(crate) const AU_KM: f64 = 149_597_870.7;
94
95/// Illumination geometry of the Moon at a Julian Ephemeris Day (TT), from the
96/// Moon's and Sun's geocentric positions (ELP/MPP02 + VSOP87D). Computed per
97/// Meeus, *Astronomical Algorithms* 2e, chapter 48 — using the true geocentric
98/// elongation (including the Moon's ecliptic latitude) and the Sun/Moon
99/// distances, not a longitude-only approximation.
100#[derive(Debug, Clone, Copy, PartialEq)]
101pub struct MoonPhase {
102    /// Geocentric ecliptic-longitude elongation Moon−Sun, degrees `[0, 360)`:
103    /// 0 = new, 90 = first quarter, 180 = full, 270 = last quarter. Drives the
104    /// phase name and the waxing/waning sense.
105    pub elongation_deg: f64,
106    /// Phase angle *i* (the Sun–Moon–Earth angle), degrees `[0, 180]`: 0 = full,
107    /// 180 = new (Meeus eq. 48.3).
108    pub phase_angle_deg: f64,
109    /// Illuminated fraction of the Moon's disc, `0.0`–`1.0` (Meeus eq. 48.1,
110    /// `k = (1 + cos i) / 2`).
111    pub illuminated_fraction: f64,
112    /// `true` while waxing (new → full, elongation in `(0, 180)`), `false` while
113    /// waning (full → new).
114    pub waxing: bool,
115}
116
117/// Compute the Moon's phase geometry (elongation, phase angle, illuminated
118/// fraction) at the given Julian Ephemeris Day in Terrestrial Time.
119///
120/// Uses the true geocentric elongation ψ from Meeus eq. 48.2
121/// (`cos ψ = cos β · cos(λ_moon − λ_sun)`, β = the Moon's ecliptic latitude) and
122/// the phase angle from eq. 48.3 with the real Sun and Moon distances, so the
123/// illuminated fraction is accurate near the quarters where the latitude term
124/// matters — not the longitude-only shortcut.
125#[must_use]
126pub fn moon_phase(jde_tt: f64) -> MoonPhase {
127    let moon = moon_position(jde_tt);
128    let sun = solar_ecliptic_state(jde_tt);
129
130    let lam_m = moon.longitude_degrees.to_radians();
131    let lam_s = sun.apparent_longitude_degrees.to_radians();
132    let beta = moon.latitude_degrees.to_radians();
133
134    // Meeus 48.2: true geocentric elongation ψ (the Sun's ecliptic latitude ≈ 0).
135    let cos_psi = beta.cos() * (lam_m - lam_s).cos();
136    let psi = cos_psi.clamp(-1.0, 1.0).acos(); // [0, π]
137
138    // Meeus 48.3: phase angle i from Sun distance R and Moon distance Δ.
139    let r_km = sun.radius_au * AU_KM;
140    let delta_km = moon.distance_km;
141    let i = (r_km * psi.sin()).atan2(delta_km - r_km * cos_psi); // [0, π]
142    let illuminated_fraction = (1.0 + i.cos()) / 2.0;
143
144    // Longitude-difference elongation for phase ordering / waxing sense.
145    let elongation_deg =
146        (moon.longitude_degrees - sun.apparent_longitude_degrees).rem_euclid(360.0);
147
148    MoonPhase {
149        elongation_deg,
150        phase_angle_deg: i.to_degrees(),
151        illuminated_fraction,
152        waxing: elongation_deg < 180.0,
153    }
154}
155
156/// Evaluate a VSOP87 series: a sum over powers of `tau`, each power weighting a
157/// sum of `A * cos(B + C * tau)` terms.
158fn eval_vsop_series(series: &[&[[f64; 3]]], tau: f64) -> f64 {
159    let mut result = 0.0;
160    let mut tau_pow = 1.0;
161    for terms in series {
162        let mut sum = 0.0;
163        for term in *terms {
164            sum += term[0] * (term[1] + term[2] * tau).cos();
165        }
166        result += sum * tau_pow;
167        tau_pow *= tau;
168    }
169    result
170}
171
172/// The five Delaunay fundamental arguments (radians) at Julian century `t` (TT).
173/// Returns `(l, l', F, D, Ω)`. Source: IERS Conventions (2010), Table 5.2a.
174pub(crate) fn delaunay_args(t: f64) -> (f64, f64, f64, f64, f64) {
175    let t2 = t * t;
176    let t3 = t2 * t;
177    let t4 = t3 * t;
178    // `% 1296000.0` (arcseconds in 360°) reduces each argument; the f64 `%`
179    // operator matches JavaScript's truncated remainder for a faithful port.
180    let l = ((485868.249036 + 1717915923.2178 * t + 31.8792 * t2 + 0.051635 * t3
181        - 0.00024470 * t4)
182        % 1296000.0)
183        * ARCSEC_TO_RAD;
184    let lp = ((1287104.79305 + 129596581.0481 * t - 0.5532 * t2 + 0.000136 * t3 - 0.00001149 * t4)
185        % 1296000.0)
186        * ARCSEC_TO_RAD;
187    let f = ((335779.526232 + 1739527262.8478 * t - 12.7512 * t2 - 0.001037 * t3
188        + 0.00000417 * t4)
189        % 1296000.0)
190        * ARCSEC_TO_RAD;
191    let d = ((1072260.70369 + 1602961601.2090 * t - 6.3706 * t2 + 0.006593 * t3 - 0.00003169 * t4)
192        % 1296000.0)
193        * ARCSEC_TO_RAD;
194    let om = ((450160.398036 - 6962890.5431 * t + 7.4722 * t2 + 0.007702 * t3 - 0.00005939 * t4)
195        % 1296000.0)
196        * ARCSEC_TO_RAD;
197    (l, lp, f, d, om)
198}
199
200/// Nutation in longitude (Δψ) in arcseconds, IAU2000B (77 lunisolar terms).
201pub(crate) fn nutation_dpsi(l: f64, lp: f64, f: f64, d: f64, om: f64, t: f64) -> f64 {
202    let mut dpsi = 0.0;
203    for row in NUT_COEFFS {
204        let arg = row[0] * l + row[1] * lp + row[2] * f + row[3] * d + row[4] * om;
205        dpsi += (row[5] + row[6] * t) * arg.sin();
206    }
207    // 0.1 microarcseconds -> arcseconds.
208    dpsi / 1e7
209}
210
211/// Nutation in obliquity (Δε) in arcseconds, IAU2000B (parallel to `NUT_COEFFS`).
212pub(crate) fn nutation_deps(l: f64, lp: f64, f: f64, d: f64, om: f64, t: f64) -> f64 {
213    let mut deps = 0.0;
214    for (row, obliq) in NUT_COEFFS.iter().zip(NUT_OBLIQ) {
215        let arg = row[0] * l + row[1] * lp + row[2] * f + row[3] * d + row[4] * om;
216        deps += (obliq[0] + obliq[1] * t) * arg.cos();
217    }
218    // 0.1 microarcseconds -> arcseconds.
219    deps / 1e7
220}
221
222/// Normalize an angle in radians to `[0, 2π)`.
223pub(crate) fn normalize_radians(rad: f64) -> f64 {
224    ((rad % TAU) + TAU) % TAU
225}
stem_branch/lib.rs

stem_branch/
lib.rs
