sidereon_core/astro/almanac/
mod.rs

1//! Astronomical almanac event finders built on the shared event finder.
2
3use std::cell::RefCell;
4use std::fmt;
5
6use crate::astro::apparent::apparent_geocentric;
7use crate::astro::events::{EventFinder, EventFinderError};
8use crate::astro::frames::transforms::{geodetic_to_itrs, FrameTransformError, GeodeticStationKm};
9use crate::astro::math::vec3::{dot3, norm3};
10use crate::astro::passes::UtcInstant;
11use crate::astro::spk::{Spk, SpkError};
12use crate::astro::{
13    constants::{
14        time::{SECONDS_PER_DAY, SECONDS_PER_HOUR},
15        units::MICROSECONDS_PER_SECOND,
16    },
17    events::CrossingEvent,
18};
19use crate::validate;
20
21mod eclipse;
22mod ecliptic;
23mod phases;
24mod planets;
25mod seasons;
26#[cfg(test)]
27mod tests;
28mod transits;
29
30pub use eclipse::lunar_solar_eclipses;
31pub use ecliptic::{geocentric_ecliptic, EclipticLonLat};
32pub use phases::{moon_phase_deg, moon_phases};
33pub use planets::planetary_events;
34pub use seasons::seasons;
35pub use transits::meridian_transits;
36
37pub(crate) const NAIF_SUN: i32 = 10;
38pub(crate) const NAIF_MOON: i32 = 301;
39
40pub(crate) const SEASON_PLANET_STEP_MAX_SECONDS: f64 = SECONDS_PER_DAY;
41pub(crate) const PHASE_STEP_MAX_SECONDS: f64 = 3.0 * SECONDS_PER_DAY;
42pub(crate) const TRANSIT_STEP_MAX_SECONDS: f64 = SECONDS_PER_HOUR;
43
44/// Which ephemeris backs an almanac computation.
45#[derive(Clone, Copy)]
46pub enum EphemerisSource<'a> {
47    /// DE-series kernel already loaded by the caller.
48    Spk(&'a Spk),
49    /// Analytic Sun/Moon series.
50    Analytic,
51}
52
53/// Seasonal marker.
54#[non_exhaustive]
55#[derive(Debug, Clone, Copy, PartialEq, Eq)]
56pub enum SeasonKind {
57    MarchEquinox,
58    JuneSolstice,
59    SeptemberEquinox,
60    DecemberSolstice,
61}
62
63/// Principal lunar phase.
64#[non_exhaustive]
65#[derive(Debug, Clone, Copy, PartialEq, Eq)]
66pub enum MoonPhaseKind {
67    New,
68    FirstQuarter,
69    Full,
70    LastQuarter,
71}
72
73/// Planetary ecliptic-longitude event.
74#[non_exhaustive]
75#[derive(Debug, Clone, Copy, PartialEq, Eq)]
76pub enum PlanetaryEventKind {
77    Conjunction,
78    Opposition,
79}
80
81/// Meridian transit kind.
82#[non_exhaustive]
83#[derive(Debug, Clone, Copy, PartialEq, Eq)]
84pub enum CulminationKind {
85    Upper,
86    Lower,
87}
88
89/// Lunar or solar eclipse kind.
90#[non_exhaustive]
91#[derive(Debug, Clone, Copy, PartialEq, Eq)]
92pub enum EclipseKind {
93    LunarPenumbral,
94    LunarPartial,
95    LunarTotal,
96    SolarPartial,
97    SolarAnnular,
98    SolarTotal,
99    SolarHybrid,
100}
101
102/// Planet selector for almanac event finders.
103#[non_exhaustive]
104#[derive(Debug, Clone, Copy, PartialEq, Eq)]
105pub enum Planet {
106    Mercury,
107    Venus,
108    Mars,
109    Jupiter,
110    Saturn,
111    Uranus,
112    Neptune,
113}
114
115/// Body selector for meridian transits.
116#[non_exhaustive]
117#[derive(Debug, Clone, Copy, PartialEq, Eq)]
118pub enum TransitBody {
119    Sun,
120    Moon,
121    Planet(Planet),
122}
123
124/// One seasonal marker event.
125#[non_exhaustive]
126#[derive(Debug, Clone, Copy, PartialEq, Eq)]
127pub struct SeasonEvent {
128    pub time: UtcInstant,
129    pub kind: SeasonKind,
130}
131
132/// One lunar phase event.
133#[non_exhaustive]
134#[derive(Debug, Clone, Copy, PartialEq, Eq)]
135pub struct MoonPhaseEvent {
136    pub time: UtcInstant,
137    pub kind: MoonPhaseKind,
138}
139
140/// One planetary opposition or conjunction.
141#[non_exhaustive]
142#[derive(Debug, Clone, Copy, PartialEq)]
143pub struct PlanetaryEvent {
144    pub time: UtcInstant,
145    pub planet: Planet,
146    pub kind: PlanetaryEventKind,
147    pub elongation_deg: f64,
148}
149
150/// One upper or lower culmination.
151#[non_exhaustive]
152#[derive(Debug, Clone, Copy, PartialEq)]
153pub struct CulminationEvent {
154    pub time: UtcInstant,
155    pub kind: CulminationKind,
156    pub altitude_deg: f64,
157}
158
159/// One lunar or solar eclipse event.
160#[non_exhaustive]
161#[derive(Debug, Clone, Copy, PartialEq)]
162pub struct EclipseEvent {
163    pub time_maximum: UtcInstant,
164    pub kind: EclipseKind,
165    pub magnitude: f64,
166    pub moon_latitude_deg: f64,
167    pub gamma: f64,
168    pub uncertain: bool,
169}
170
171/// Error returned by almanac computations.
172#[non_exhaustive]
173#[derive(Debug, Clone, PartialEq)]
174pub enum AlmanacError {
175    Finder(EventFinderError),
176    Spk(SpkError),
177    Frame(&'static str),
178    EphemerisRequired,
179    InferiorPlanetOpposition,
180    InvalidInput {
181        field: &'static str,
182        reason: &'static str,
183    },
184}
185
186impl fmt::Display for AlmanacError {
187    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
188        match self {
189            Self::Finder(error) => write!(f, "event finder failed: {error}"),
190            Self::Spk(error) => write!(f, "SPK failed: {error}"),
191            Self::Frame(label) => write!(f, "frame reduction failed: {label}"),
192            Self::EphemerisRequired => write!(f, "SPK ephemeris is required"),
193            Self::InferiorPlanetOpposition => {
194                write!(f, "opposition is not defined for an inferior planet")
195            }
196            Self::InvalidInput { field, reason } => {
197                write!(f, "invalid almanac input {field}: {reason}")
198            }
199        }
200    }
201}
202
203impl std::error::Error for AlmanacError {}
204
205pub(crate) fn validate_scan_controls(
206    step_seconds: f64,
207    time_tolerance_seconds: f64,
208    step_max_seconds: f64,
209) -> Result<(), AlmanacError> {
210    validate::positive_step(step_seconds, "step_seconds").map_err(map_field_error)?;
211    validate::positive_step(time_tolerance_seconds, "time_tolerance_seconds")
212        .map_err(map_field_error)?;
213    if step_seconds > step_max_seconds {
214        return Err(AlmanacError::InvalidInput {
215            field: "step_seconds",
216            reason: "exceeds maximum",
217        });
218    }
219    Ok(())
220}
221
222pub(crate) fn validate_station(station: &GeodeticStationKm) -> Result<(), AlmanacError> {
223    geodetic_to_itrs(
224        station.latitude_deg,
225        station.longitude_deg,
226        station.altitude_km,
227    )
228    .map(|_| ())
229    .map_err(map_frame_input)
230}
231
232pub(crate) fn event_finder(
233    start: UtcInstant,
234    end: UtcInstant,
235    step_seconds: f64,
236    time_tolerance_seconds: f64,
237) -> Result<EventFinder, AlmanacError> {
238    EventFinder::new(
239        0.0,
240        seconds_between(start, end)?,
241        step_seconds,
242        time_tolerance_seconds,
243    )
244    .map_err(AlmanacError::Finder)
245}
246
247pub(crate) fn seconds_between(start: UtcInstant, end: UtcInstant) -> Result<f64, AlmanacError> {
248    let span = end
249        .unix_microseconds()
250        .checked_sub(start.unix_microseconds())
251        .ok_or(AlmanacError::Finder(EventFinderError::InvalidInput {
252            field: "time_window",
253            reason: "start/end span overflows i64 microseconds",
254        }))?;
255    Ok(span as f64 / MICROSECONDS_PER_SECOND)
256}
257
258pub(crate) fn instant_at_offset_seconds(start: UtcInstant, offset_seconds: f64) -> UtcInstant {
259    UtcInstant::from_unix_microseconds(
260        start.unix_microseconds() + (offset_seconds * MICROSECONDS_PER_SECOND).floor() as i64,
261    )
262}
263
264pub(crate) fn offset_instant(start: UtcInstant, offset_seconds: f64) -> UtcInstant {
265    instant_at_offset_seconds(start, offset_seconds)
266}
267
268pub(crate) fn body_ecliptic(
269    source: EphemerisSource<'_>,
270    target_naif: i32,
271    time: UtcInstant,
272) -> Result<EclipticLonLat, AlmanacError> {
273    let ts = time.time_scales();
274    let pos = apparent_geocentric(target_naif, &ts, source)?;
275    geocentric_ecliptic(pos, &ts)
276}
277
278pub(crate) fn apparent_km(
279    source: EphemerisSource<'_>,
280    target_naif: i32,
281    time: UtcInstant,
282) -> Result<[f64; 3], AlmanacError> {
283    let pos_m = apparent_geocentric(target_naif, &time.time_scales(), source)?;
284    Ok([pos_m[0] * 1.0e-3, pos_m[1] * 1.0e-3, pos_m[2] * 1.0e-3])
285}
286
287pub(crate) fn find_angle_crossing_times<F>(
288    start: UtcInstant,
289    end: UtcInstant,
290    step_seconds: f64,
291    time_tolerance_seconds: f64,
292    target_deg: f64,
293    angle_fn: F,
294) -> Result<Vec<UtcInstant>, AlmanacError>
295where
296    F: Fn(UtcInstant) -> Result<f64, AlmanacError>,
297{
298    let finder = event_finder(start, end, step_seconds, time_tolerance_seconds)?;
299    let latch = RefCell::new(None);
300    let crossings = finder
301        .find_crossings(
302            |offset_seconds| {
303                latch_scalar(&latch, || {
304                    let time = instant_at_offset_seconds(start, offset_seconds);
305                    let angle = angle_fn(time)?;
306                    Ok((angle - target_deg).to_radians().sin())
307                })
308            },
309            0.0,
310        )
311        .map_err(|error| latched_or_finder(error, &latch))?;
312
313    let mut times = Vec::new();
314    for crossing in crossings {
315        let time = instant_at_offset_seconds(start, crossing.time_seconds);
316        let angle = angle_fn(time)?;
317        if (angle - target_deg).to_radians().cos() > 0.0 {
318            times.push(time);
319        }
320    }
321    Ok(times)
322}
323
324pub(crate) fn latch_scalar<F>(latch: &RefCell<Option<AlmanacError>>, f: F) -> f64
325where
326    F: FnOnce() -> Result<f64, AlmanacError>,
327{
328    match f() {
329        Ok(value) if value.is_finite() => value,
330        Ok(_) => {
331            latch_error(
332                latch,
333                AlmanacError::InvalidInput {
334                    field: "predicate",
335                    reason: "not finite",
336                },
337            );
338            f64::NAN
339        }
340        Err(error) => {
341            latch_error(latch, error);
342            f64::NAN
343        }
344    }
345}
346
347pub(crate) fn latched_or_finder(
348    error: EventFinderError,
349    latch: &RefCell<Option<AlmanacError>>,
350) -> AlmanacError {
351    latch
352        .borrow()
353        .clone()
354        .unwrap_or(AlmanacError::Finder(error))
355}
356
357pub(crate) fn latch_error(latch: &RefCell<Option<AlmanacError>>, error: AlmanacError) {
358    if latch.borrow().is_none() {
359        *latch.borrow_mut() = Some(error);
360    }
361}
362
363pub(crate) fn crossing_time(start: UtcInstant, crossing: CrossingEvent) -> UtcInstant {
364    instant_at_offset_seconds(start, crossing.time_seconds)
365}
366
367pub(crate) fn planet_naif(planet: Planet) -> i32 {
368    match planet {
369        Planet::Mercury => 1,
370        Planet::Venus => 2,
371        Planet::Mars => 4,
372        Planet::Jupiter => 5,
373        Planet::Saturn => 6,
374        Planet::Uranus => 7,
375        Planet::Neptune => 8,
376    }
377}
378
379pub(crate) fn transit_body_naif(body: TransitBody) -> i32 {
380    match body {
381        TransitBody::Sun => NAIF_SUN,
382        TransitBody::Moon => NAIF_MOON,
383        TransitBody::Planet(planet) => planet_naif(planet),
384    }
385}
386
387pub(crate) fn is_inferior(planet: Planet) -> bool {
388    matches!(planet, Planet::Mercury | Planet::Venus)
389}
390
391pub(crate) fn wrap360(degrees: f64) -> f64 {
392    degrees.rem_euclid(360.0)
393}
394
395pub(crate) fn angular_separation_rad(a: [f64; 3], b: [f64; 3]) -> Result<f64, AlmanacError> {
396    let na = norm_checked(a, "a")?;
397    let nb = norm_checked(b, "b")?;
398    let cos_sep = (dot3(a, b) / (na * nb)).clamp(-1.0, 1.0);
399    Ok(cos_sep.acos())
400}
401
402pub(crate) fn norm_checked(vector: [f64; 3], field: &'static str) -> Result<f64, AlmanacError> {
403    if vector.iter().any(|value| !value.is_finite()) {
404        return Err(AlmanacError::InvalidInput {
405            field,
406            reason: "components must be finite",
407        });
408    }
409    let norm = norm3(vector);
410    if !norm.is_finite() {
411        return Err(AlmanacError::InvalidInput {
412            field,
413            reason: "norm must be finite",
414        });
415    }
416    if norm == 0.0 {
417        return Err(AlmanacError::InvalidInput {
418            field,
419            reason: "degenerate",
420        });
421    }
422    Ok(norm)
423}
424
425pub(crate) fn map_field_error(error: validate::FieldError) -> AlmanacError {
426    AlmanacError::InvalidInput {
427        field: error.field(),
428        reason: error.reason(),
429    }
430}
431
432fn map_frame_input(error: FrameTransformError) -> AlmanacError {
433    let FrameTransformError::InvalidInput { field, reason } = error;
434    AlmanacError::InvalidInput { field, reason }
435}
sidereon_core/astro/almanac/mod.rs

sidereon_core/astro/almanac/
mod.rs
