pleiades_eclipse/

engine.rs

//! The public eclipse search engine.

use crate::ephemeris::{apparent_sun_longitude_deg, sample_sun_moon};
use crate::error::{EclipseError, WINDOW_END_JD, WINDOW_START_JD};
use crate::geometry::{classify_lunar, classify_solar, sub_shadow_point};
use crate::local::{is_locally_visible, local_circumstances_for, LocalCircumstances};
use crate::saros::saros_series;
use crate::syzygy::{find_syzygies, Syzygy, STEP_DAYS};
use crate::types::{Eclipse, EclipseFilter, EclipseKind, EclipseType, Node};
use pleiades_apparent::Atmosphere;
use pleiades_backend::EphemerisBackend;
use pleiades_types::{Instant, JulianDay, Longitude, ObserverLocation, TimeScale};

/// Backstop cap on how many global eclipses `next/previous_local_eclipse`
/// inspect before returning `None` (a locally-visible eclipse always occurs
/// far sooner; this only guards against a pathological non-terminating walk).
const MAX_LOCAL_SEARCH: usize = 4000;

/// Half-width of the bracket around a syzygy in which `refine_greatest` looks
/// for greatest eclipse. `eclipses_in_range` widens its syzygy scan by the same
/// amount so that an eclipse whose greatest instant lies inside the requested
/// range is found even when its syzygy falls just outside it.
const GREATEST_BRACKET_DAYS: f64 = 0.25;

/// Width of the first span `next_eclipse`/`previous_eclipse` search outward
/// from the query instant; later spans double. Under every filter the admitted
/// eclipses recur well within a year, so away from the window edges the first
/// span holds the answer; the value only tunes cost, never the result.
const FIRST_SEARCH_SPAN_DAYS: f64 = 400.0;

/// Searches for global/geocentric eclipses over a chosen [`EphemerisBackend`].
///
/// Backed by the packaged data, results are valid only within the
/// 1900-01-01..2100-01-01 window (see [`crate`]); out-of-window requests fail
/// closed with [`EclipseError::OutOfWindow`].
pub struct EclipseEngine<B> {
    backend: B,
}

impl<B: EphemerisBackend> EclipseEngine<B> {
    /// Creates an engine that draws Sun/Moon positions from `backend`.
    pub fn new(backend: B) -> Self {
        Self { backend }
    }

    /// Returns every eclipse admitted by `filter` with greatest eclipse in
    /// `[start, end]` (both bounds inclusive), in chronological order. Fails
    /// closed if either bound is outside the supported window.
    ///
    /// Selection is by the greatest-eclipse instant, not by the new/full moon
    /// that precedes or follows it by minutes, and that instant does not depend
    /// on the range searched. Adjacent ranges that meet at an instant therefore
    /// partition the eclipses between them, except that an eclipse whose
    /// greatest instant falls exactly on the shared bound appears in both.
    pub fn eclipses_in_range(
        &self,
        start: Instant,
        end: Instant,
        filter: EclipseFilter,
    ) -> Result<Vec<Eclipse>, EclipseError> {
        let start_jd = start.julian_day.days();
        let end_jd = end.julian_day.days();
        self.check_window(start_jd)?;
        self.check_window(end_jd)?;

        // Greatest eclipse lies within GREATEST_BRACKET_DAYS of its syzygy, so a
        // syzygy up to that far outside `[start, end]` can still carry an eclipse
        // whose greatest instant is inside it. Scan that wider span, then select
        // on the greatest instant below.
        //
        // The syzygy scanner (`find_one`) samples a grid of STEP_DAYS multiples
        // that may begin one step before its requested start and probes one step
        // past its requested end, and `sample_sun_moon` issues a light-time-
        // retarded Sun query ~0.006 days before the nominal epoch. Together these
        // mean:
        //   - At the END: the scanner queries up to `scan_end + STEP_DAYS`; the
        //     packaged backend has no data beyond WINDOW_END_JD, so we clamp
        //     `scan_end` to `WINDOW_END_JD - STEP_DAYS`.
        //   - At the START: WINDOW_START_JD is itself a multiple of STEP_DAYS, so
        //     clamping `scan_start` to `WINDOW_START_JD + STEP_DAYS` keeps the
        //     first grid point at or after it, and the retarded query (~0.006 d
        //     earlier) within coverage.
        // Both clamps are safe: no corpus eclipse falls within 0.5 d of either bound.
        let scan_start = (start_jd - GREATEST_BRACKET_DAYS).max(WINDOW_START_JD + STEP_DAYS);
        let scan_end = (end_jd + GREATEST_BRACKET_DAYS).min(WINDOW_END_JD - STEP_DAYS);

        let mut out = Vec::new();
        for event in find_syzygies(&self.backend, scan_start, scan_end)? {
            let Some(eclipse) = self.build(event.syzygy, event.julian_day)? else {
                continue;
            };
            let greatest_jd = eclipse.greatest_eclipse.julian_day.days();
            if filter.admits(eclipse.kind) && (start_jd..=end_jd).contains(&greatest_jd) {
                out.push(eclipse);
            }
        }
        Ok(out)
    }

    /// Returns the first eclipse admitted by `filter` whose greatest eclipse is
    /// strictly after `after`, or `None` if none remains before the window end.
    /// Fails closed if `after` is outside the supported window.
    ///
    /// The cost grows with the distance to the eclipse found, not to the window
    /// end.
    pub fn next_eclipse(
        &self,
        after: Instant,
        filter: EclipseFilter,
    ) -> Result<Option<Eclipse>, EclipseError> {
        self.search_forward(after.julian_day.days(), FIRST_SEARCH_SPAN_DAYS, filter)
    }

    /// Returns the last eclipse admitted by `filter` whose greatest eclipse is
    /// strictly before `before`, or `None` if none exists after the window start.
    /// Fails closed if `before` is outside the supported window.
    ///
    /// The cost grows with the distance to the eclipse found, not to the window
    /// start.
    pub fn previous_eclipse(
        &self,
        before: Instant,
        filter: EclipseFilter,
    ) -> Result<Option<Eclipse>, EclipseError> {
        self.search_backward(before.julian_day.days(), FIRST_SEARCH_SPAN_DAYS, filter)
    }

    // Outward search behind `next_eclipse`/`previous_eclipse`.
    //
    // Scans adjacent spans moving away from the query instant (widths
    // `first_span`, 2×, 4×, …, the last clamped to the window edge) and returns
    // the first admitted eclipse past the query instant in the first span that
    // has one. This is exactly the eclipse a single scan from the query instant
    // to the window edge returns first: `eclipses_in_range` selects by the
    // greatest-eclipse instant inside its inclusive bounds, and that instant does
    // not depend on the range searched (the syzygy grid is anchored to absolute
    // STEP_DAYS multiples), so the spans together cover the same instants with
    // the same eclipses and the earlier, empty spans hold nothing nearer. An
    // eclipse exactly on a shared span bound shows up in both spans; harmless,
    // since the strict comparison against the query instant decides and the
    // nearer span is searched first.

    fn search_forward(
        &self,
        after_jd: f64,
        first_span: f64,
        filter: EclipseFilter,
    ) -> Result<Option<Eclipse>, EclipseError> {
        self.check_window(after_jd)?;
        let mut near = after_jd;
        let mut span = first_span;
        loop {
            // `eclipses_in_range` clamps its scan end to WINDOW_END_JD - STEP_DAYS,
            // so ending the last span at WINDOW_END_JD itself is safe and correct.
            let far = (near + span).min(WINDOW_END_JD);
            let found = self
                .eclipses_in_range(tdb(near), tdb(far), filter)?
                .into_iter()
                .find(|e| e.greatest_eclipse.julian_day.days() > after_jd);
            // `span` stays positive and grows, so `far` strictly advances until
            // it reaches the window end, which ends the loop.
            if found.is_some() || far >= WINDOW_END_JD {
                return Ok(found);
            }
            near = far;
            span *= 2.0;
        }
    }

    fn search_backward(
        &self,
        before_jd: f64,
        first_span: f64,
        filter: EclipseFilter,
    ) -> Result<Option<Eclipse>, EclipseError> {
        self.check_window(before_jd)?;
        let mut near = before_jd;
        let mut span = first_span;
        loop {
            let far = (near - span).max(WINDOW_START_JD);
            let found = self
                .eclipses_in_range(tdb(far), tdb(near), filter)?
                .into_iter()
                .rev()
                .find(|e| e.greatest_eclipse.julian_day.days() < before_jd);
            if found.is_some() || far <= WINDOW_START_JD {
                return Ok(found);
            }
            near = far;
            span *= 2.0;
        }
    }

    /// Local (per-observer) circumstances for an already-found `eclipse`.
    ///
    /// Returns full circumstances even when the eclipse is not visible from
    /// `observer` (all contacts below the horizon); inspect `any_phase_visible`
    /// (via the returned variant) to test visibility. Solar contact instants are
    /// observer-dependent (topocentric); lunar contact instants are global with
    /// per-observer visibility.
    ///
    /// `atmosphere` is `pleiades_apparent::Atmosphere` from the
    /// `pleiades-apparent` release this crate pins (0.7 and later); build it
    /// from that same release.
    pub fn local_circumstances(
        &self,
        eclipse: &Eclipse,
        observer: &ObserverLocation,
        atmosphere: Atmosphere,
    ) -> Result<LocalCircumstances, EclipseError> {
        observer
            .validate()
            .map_err(|e| EclipseError::InvalidObserver {
                detail: e.to_string(),
            })?;
        check_atmosphere(atmosphere)?;
        local_circumstances_for(&self.backend, eclipse, observer, atmosphere)
    }

    /// The next eclipse admitted by `filter`, strictly after `after`, that is
    /// locally visible from `observer` (any phase above the horizon), paired with
    /// its local circumstances. Walks the global `next_eclipse` sequence and
    /// returns the first locally-visible one, so the result is a strict refinement
    /// of the global engine.
    pub fn next_local_eclipse(
        &self,
        after: Instant,
        observer: &ObserverLocation,
        filter: EclipseFilter,
        atmosphere: Atmosphere,
    ) -> Result<Option<(Eclipse, LocalCircumstances)>, EclipseError> {
        observer
            .validate()
            .map_err(|e| EclipseError::InvalidObserver {
                detail: e.to_string(),
            })?;
        check_atmosphere(atmosphere)?;
        let mut cursor = after;
        // Bounded walk: no more than MAX_LOCAL_SEARCH global eclipses inspected
        // before giving up (backstop; a locally-visible eclipse always occurs well
        // within the window). ~2 eclipses/year × 200 yr ≈ 1200 global eclipses max.
        for _ in 0..MAX_LOCAL_SEARCH {
            let Some(eclipse) = self.next_eclipse(cursor, filter)? else {
                return Ok(None);
            };
            let local = local_circumstances_for(&self.backend, &eclipse, observer, atmosphere)?;
            if is_locally_visible(&local) {
                return Ok(Some((eclipse, local)));
            }
            cursor = eclipse.greatest_eclipse;
        }
        Ok(None)
    }

    /// The previous eclipse admitted by `filter`, strictly before `before`, that
    /// is locally visible from `observer`, paired with its local circumstances.
    pub fn previous_local_eclipse(
        &self,
        before: Instant,
        observer: &ObserverLocation,
        filter: EclipseFilter,
        atmosphere: Atmosphere,
    ) -> Result<Option<(Eclipse, LocalCircumstances)>, EclipseError> {
        observer
            .validate()
            .map_err(|e| EclipseError::InvalidObserver {
                detail: e.to_string(),
            })?;
        check_atmosphere(atmosphere)?;
        let mut cursor = before;
        for _ in 0..MAX_LOCAL_SEARCH {
            let Some(eclipse) = self.previous_eclipse(cursor, filter)? else {
                return Ok(None);
            };
            let local = local_circumstances_for(&self.backend, &eclipse, observer, atmosphere)?;
            if is_locally_visible(&local) {
                return Ok(Some((eclipse, local)));
            }
            cursor = eclipse.greatest_eclipse;
        }
        Ok(None)
    }

    fn check_window(&self, jd: f64) -> Result<(), EclipseError> {
        if !(WINDOW_START_JD..=WINDOW_END_JD).contains(&jd) {
            Err(EclipseError::OutOfWindow { julian_day: jd })
        } else {
            Ok(())
        }
    }

    fn build(&self, syzygy: Syzygy, syzygy_jd: f64) -> Result<Option<Eclipse>, EclipseError> {
        let greatest_jd = self.refine_greatest(syzygy, syzygy_jd)?;
        let sample = sample_sun_moon(&self.backend, greatest_jd)?;
        let greatest_eclipse = Instant::new(JulianDay::from_days(greatest_jd), TimeScale::Tdb);
        // Compute the apparent geocentric solar longitude of date for eclipsed_longitude.
        // Geometric sampling (separation, classification) uses the Mean frame — that is
        // correct and unchanged. But eclipsed_longitude must be apparent-of-date per the
        // spec (Task 10 gate: ≤1 arcsecond); mean would be ~20–25″ off.
        let apparent_sun_lon = apparent_sun_longitude_deg(&self.backend, greatest_jd)?;
        let eclipsed_longitude = match syzygy {
            Syzygy::NewMoon => Longitude::from_degrees(apparent_sun_lon),
            // For lunar eclipses the eclipsed body is the Moon, which is opposite the Sun;
            // eclipsed_longitude is the apparent solar longitude + 180° (corpus MANIFEST).
            Syzygy::FullMoon => Longitude::from_degrees(apparent_sun_lon + 180.0),
        };
        // Node: ascending (North) if the Moon's latitude is increasing through 0.
        let later = sample_sun_moon(&self.backend, greatest_jd + 0.01)?;
        let near_node = if later.moon_latitude_deg >= sample.moon_latitude_deg {
            Node::North
        } else {
            Node::South
        };

        let eclipse = match syzygy {
            Syzygy::NewMoon => {
                let Some(c) = classify_solar(&sample) else {
                    return Ok(None);
                };
                Eclipse {
                    kind: EclipseKind::Solar,
                    eclipse_type: EclipseType::Solar(c.eclipse_type),
                    greatest_eclipse,
                    magnitude: c.magnitude,
                    gamma: c.gamma,
                    saros_series: saros_series(EclipseKind::Solar, greatest_jd),
                    eclipsed_longitude,
                    near_node,
                    greatest_eclipse_location: Some(sub_shadow_point(&sample, greatest_jd)),
                }
            }
            Syzygy::FullMoon => {
                let Some(c) = classify_lunar(&sample) else {
                    return Ok(None);
                };
                Eclipse {
                    kind: EclipseKind::Lunar,
                    eclipse_type: EclipseType::Lunar(c.eclipse_type),
                    greatest_eclipse,
                    magnitude: c.magnitude,
                    gamma: c.gamma,
                    saros_series: saros_series(EclipseKind::Lunar, greatest_jd),
                    eclipsed_longitude,
                    near_node,
                    greatest_eclipse_location: None,
                }
            }
        };
        Ok(Some(eclipse))
    }

    /// Golden-section minimize the Sun–Moon (or Moon–antisolar) separation in a
    /// ±GREATEST_BRACKET_DAYS bracket around the syzygy to find greatest eclipse.
    fn refine_greatest(&self, syzygy: Syzygy, syzygy_jd: f64) -> Result<f64, EclipseError> {
        use crate::geometry::separation_for;
        let phi = 0.618_033_988_75_f64;
        let (mut a, mut b) = (
            syzygy_jd - GREATEST_BRACKET_DAYS,
            syzygy_jd + GREATEST_BRACKET_DAYS,
        );
        let mut c = b - (b - a) * phi;
        let mut d = a + (b - a) * phi;
        let mut fc = separation_for(syzygy, &sample_sun_moon(&self.backend, c)?);
        let mut fd = separation_for(syzygy, &sample_sun_moon(&self.backend, d)?);
        while (b - a) > 0.5 / 86_400.0 {
            if fc < fd {
                b = d;
                d = c;
                fd = fc;
                c = b - (b - a) * phi;
                fc = separation_for(syzygy, &sample_sun_moon(&self.backend, c)?);
            } else {
                a = c;
                c = d;
                fc = fd;
                d = a + (b - a) * phi;
                fd = separation_for(syzygy, &sample_sun_moon(&self.backend, d)?);
            }
        }
        Ok(0.5 * (a + b))
    }
}

fn tdb(jd: f64) -> Instant {
    Instant::new(JulianDay::from_days(jd), TimeScale::Tdb)
}

fn check_atmosphere(atmos: Atmosphere) -> Result<(), EclipseError> {
    if !atmos.pressure_mbar.is_finite() || !atmos.temperature_c.is_finite() {
        return Err(EclipseError::InvalidAtmosphere {
            detail: format!(
                "pressure={} temp={}",
                atmos.pressure_mbar, atmos.temperature_c
            ),
        });
    }
    Ok(())
}

#[cfg(test)]
mod tests {
    use super::*;
    use pleiades_backend::test_backend::LinearSunMoon;
    use pleiades_types::{Instant, JulianDay, TimeScale};

    fn at(jd: f64) -> Instant {
        Instant::new(JulianDay::from_days(jd), TimeScale::Tdb)
    }

    #[test]
    fn out_of_window_start_fails_closed() {
        let engine = EclipseEngine::new(LinearSunMoon::new_moon_at(2_451_550.0));
        let err = engine
            .eclipses_in_range(at(2_400_000.0), at(2_451_551.0), EclipseFilter::All)
            .unwrap_err();
        assert!(matches!(err, EclipseError::OutOfWindow { .. }));
    }

    #[test]
    fn filter_excludes_lunar() {
        // The on-node analytic backend yields a solar eclipse at every new moon
        // and a lunar one at every full moon; SolarOnly must drop the lunar ones.
        let engine =
            EclipseEngine::new(LinearSunMoon::new_moon_at(2_451_550.0).with_moon_latitude(0.0));
        let solar = engine
            .eclipses_in_range(at(2_451_549.0), at(2_451_551.0), EclipseFilter::SolarOnly)
            .unwrap();
        assert!(solar.iter().all(|e| e.kind == EclipseKind::Solar));
    }

    /// Real eclipses always recur within the 400-day first span, so the
    /// doubling path is exercised by starting from a 1-day span: reaching the
    /// next/previous eclipse (months away) then takes many doubled spans, and
    /// must give exactly the eclipse the default span gives in one. (The
    /// default-span results are pinned against a single `eclipses_in_range`
    /// scan in `tests/known_eclipses.rs`.)
    #[test]
    fn outward_search_result_does_not_depend_on_the_first_span() {
        let engine = EclipseEngine::new(pleiades_data::packaged_backend());
        // Near each edge: one direction walks to an eclipse months away, the
        // other runs into the clamped window edge.
        for t in [WINDOW_START_JD + 40.0, WINDOW_END_JD - 40.0] {
            for filter in [
                EclipseFilter::All,
                EclipseFilter::SolarOnly,
                EclipseFilter::LunarOnly,
            ] {
                assert_eq!(
                    engine.search_forward(t, 1.0, filter).unwrap(),
                    engine
                        .search_forward(t, FIRST_SEARCH_SPAN_DAYS, filter)
                        .unwrap(),
                    "forward from JD {t}, {filter:?}"
                );
                assert_eq!(
                    engine.search_backward(t, 1.0, filter).unwrap(),
                    engine
                        .search_backward(t, FIRST_SEARCH_SPAN_DAYS, filter)
                        .unwrap(),
                    "backward from JD {t}, {filter:?}"
                );
            }
        }
    }

    #[test]
    fn outward_search_with_no_eclipse_ends_at_the_window_edge() {
        // Off-node Moon: no syzygy is an eclipse, so the search must walk its
        // doubling spans out to the window edge and stop there with `None`.
        let engine =
            EclipseEngine::new(LinearSunMoon::new_moon_at(2_451_550.0).with_moon_latitude(5.0));
        let t = WINDOW_END_JD - 30.0;
        assert_eq!(
            engine.search_forward(t, 1.0, EclipseFilter::All).unwrap(),
            None
        );
        let t = WINDOW_START_JD + 30.0;
        assert_eq!(
            engine.search_backward(t, 1.0, EclipseFilter::All).unwrap(),
            None
        );
    }

    #[test]
    fn next_and_previous_eclipse_fail_closed_outside_the_window() {
        let engine = EclipseEngine::new(LinearSunMoon::new_moon_at(2_451_550.0));
        for jd in [WINDOW_START_JD - 1.0, WINDOW_END_JD + 1.0, f64::NAN] {
            assert!(matches!(
                engine.next_eclipse(at(jd), EclipseFilter::All),
                Err(EclipseError::OutOfWindow { .. })
            ));
            assert!(matches!(
                engine.previous_eclipse(at(jd), EclipseFilter::All),
                Err(EclipseError::OutOfWindow { .. })
            ));
        }
    }

    #[test]
    fn local_circumstances_returns_solar_for_a_solar_eclipse() {
        use pleiades_apparent::Atmosphere;
        use pleiades_types::{Latitude, Longitude, ObserverLocation};
        let engine =
            EclipseEngine::new(LinearSunMoon::new_moon_at(2_451_550.0).with_moon_latitude(0.0));
        let eclipse = engine
            .next_eclipse(at(2_451_549.0), EclipseFilter::SolarOnly)
            .unwrap()
            .expect("a solar eclipse");
        let observer = ObserverLocation::new(
            Latitude::from_degrees(0.0),
            Longitude::from_degrees(0.0),
            Some(0.0),
        );
        let local = engine
            .local_circumstances(&eclipse, &observer, Atmosphere::default())
            .unwrap();
        assert!(matches!(local, crate::LocalCircumstances::Solar(_)));
    }
}