// SPDX-License-Identifier: AGPL-3.0-only
// Copyright (C) 2026 Vallés Puig, Ramon

//! # Moon Ephemeris Cache, Chebyshev Segment Interpolation
//!
//! ## Scientific scope
//!
//! Precomputes the Moon's geocentric ecliptic Cartesian coordinates
//! (X, Y, Z) on Chebyshev nodes within fixed‑duration segments, then
//! evaluates any intermediate time via Clenshaw recurrence in
//! **O(degree)**, replacing the full ELP2000 series summation
//! (~200 µs → ~1 µs per query). A companion [`NutationCache`] stores
//! pre‑evaluated nutation triplets (Δψ, Δε, ε₀) at regular intervals and
//! linearly interpolates, removing the 77‑term IAU 2000B series from the
//! per‑query hot path. The Moon's shortest significant perturbation
//! period is ~5 days (evection); Chebyshev degree 8 over 4‑day segments
//! yields interpolation errors well below 1″ in geocentric position, far
//! smaller than the ~0.5° atmospheric refraction uncertainty at the
//! horizon.
//!
//! ## Technical scope
//!
//! Three public types:
//!
//! * [`MoonPositionCache`], Chebyshev cache for ELP2000 geocentric XYZ.
//! * [`NutationCache`], 2‑hour linear cache for IAU 2000B nutation.
//! * [`MoonAltitudeContext`], pre‑built context combining both caches
//!   plus precomputed observer ITRF position for the unified altitude
//!   pipeline at a fixed site.
//!
//! Plus the free function [`find_and_label_crossings`], a fused
//! scan + Brent + labelling pass that records crossing direction
//! directly from the bracketed sign change, eliminating the two probe
//! evaluations per crossing required by the generic `find_crossings` +
//! `label_crossings` pipeline in `intervals.rs`.
//!
//! Defaults: 4‑day segments, Chebyshev degree 8 (9 nodes), 2‑hour
//! nutation step.
//!
//! ## References
//! - Chapront‑Touzé, M. & Chapront, J. (1983). "The lunar ephemeris ELP 2000".
//!   *Astronomy & Astrophysics* 124, 50–62.
//! - Capitaine, N., Wallace, P. T., & Chapront, J. (2003). IAU 2000B
//!   nutation series. *Astronomy & Astrophysics* 412, 567–586.
//!   doi:10.1051/0004‑6361:20031539

#![allow(unreachable_pub, missing_docs)]

use crate::astro::earth_rotation::jd_ut1_from_tt_eop;
use crate::astro::earth_rotation_provider::itrs_to_equatorial_mean_j2000_rotation;
use crate::astro::nutation::nutation_iau2000b;
use crate::astro::precession::precession_matrix_iau2006;
use crate::astro::sidereal::gast_iau2006;
use crate::coordinates::centers::Geodetic;
use crate::coordinates::frames::ECEF;
use crate::coordinates::transform::context::{DefaultEop, DefaultEphemeris, DefaultNutationModel};
use crate::coordinates::transform::AstroContext;
use crate::ephemeris::Ephemeris;
use crate::qtty::*;
use crate::time::JulianDate;
use cheby;

// =============================================================================
// Constants
// =============================================================================

/// Chebyshev polynomial degree for position interpolation.
const CHEB_DEGREE: usize = 8;

/// Number of Chebyshev nodes per segment (degree + 1).
const CHEB_NODES: usize = CHEB_DEGREE + 1;

/// Duration of each Chebyshev segment in days.
const SEGMENT_DAYS: Days = Days::new(4.0);

/// J2000 mean obliquity ε₀ (IAU 2006): 84381.406″ converted to radians.
/// Used for ecliptic → equatorial rotation (constant for J2000 frame).
const J2000_OBLIQUITY_RAD: crate::qtty::Quantity<Radian> =
    crate::qtty::Quantity::<Radian>::new(84381.406 / 3600.0 * std::f64::consts::PI / 180.0);

/// Nutation cache step in days (2 hours).
const NUT_STEP_DAYS: Days = Hours::new(2.0).to_const::<Day>();

// =============================================================================
// MoonPositionCache
// =============================================================================

/// Chebyshev interpolation cache for the Moon's geocentric ecliptic
/// Cartesian coordinates (X, Y, Z) in kilometers.
///
/// The time domain is divided into segments of [`SEGMENT_DAYS`] days.
/// Within each segment the three coordinates are approximated by degree-8
/// Chebyshev polynomials fitted at the canonical nodes.
pub struct MoonPositionCache {
    /// Modified Julian Date of the first segment's start.
    mjd_start: ModifiedJulianDate,
    /// Number of segments.
    num_segments: usize,
    /// Chebyshev coefficients for X coordinate: [segment][CHEB_NODES], stored as raw km.
    cx: Vec<[f64; CHEB_NODES]>,
    /// Chebyshev coefficients for Y coordinate.
    cy: Vec<[f64; CHEB_NODES]>,
    /// Chebyshev coefficients for Z coordinate.
    cz: Vec<[f64; CHEB_NODES]>,
}

impl MoonPositionCache {
    /// Build the cache covering `[mjd_start, mjd_end]` (Modified Julian Dates).
    ///
    /// Adds a small padding on each side to accommodate Brent probes
    /// near the boundaries.
    ///
    /// # Arguments
    ///
    /// * `mjd_start`, inclusive lower bound of the cached window.
    /// * `mjd_end`, inclusive upper bound of the cached window.
    ///
    /// # Returns
    ///
    /// A populated [`MoonPositionCache`] holding Chebyshev coefficients
    /// for X, Y, Z over each segment.
    pub fn new(mjd_start: ModifiedJulianDate, mjd_end: ModifiedJulianDate) -> Self {
        let pad = Days::new(1.0); // 1-day padding on each side
        let t0 = crate::time::ModifiedJulianDate::new((mjd_start.raw() - pad).value());
        let span = mjd_end.raw() + pad - t0.raw();
        let num_segments = ((span / SEGMENT_DAYS).ceil() as usize).max(1);

        let nodes: [f64; CHEB_NODES] = cheby::nodes();
        let mut cx = Vec::with_capacity(num_segments);
        let mut cy = Vec::with_capacity(num_segments);
        let mut cz = Vec::with_capacity(num_segments);

        for seg in 0..num_segments {
            let seg_start = crate::time::ModifiedJulianDate::new(
                (t0.raw() + seg as f64 * SEGMENT_DAYS).value(),
            );
            let seg_mid = crate::time::ModifiedJulianDate::new(
                (seg_start.raw() + SEGMENT_DAYS * 0.5).value(),
            );
            let seg_half = SEGMENT_DAYS * 0.5;

            let mut vx = [0.0; CHEB_NODES];
            let mut vy = [0.0; CHEB_NODES];
            let mut vz = [0.0; CHEB_NODES];

            for k in 0..CHEB_NODES {
                let mjd_k = crate::time::ModifiedJulianDate::new(
                    (seg_mid.raw() + seg_half * nodes[k]).value(),
                );
                let pos = DefaultEphemeris::moon_geocentric(mjd_k.to::<crate::JD>());
                vx[k] = pos.x().value();
                vy[k] = pos.y().value();
                vz[k] = pos.z().value();
            }

            cx.push(cheby::fit_coeffs(&vx));
            cy.push(cheby::fit_coeffs(&vy));
            cz.push(cheby::fit_coeffs(&vz));
        }

        Self {
            mjd_start: t0,
            num_segments,
            cx,
            cy,
            cz,
        }
    }

    /// Evaluate the cached geocentric ecliptic (X, Y, Z) in km at `jd`.
    ///
    /// Falls back to full ELP2000 if `jd` is outside the cached range.
    ///
    /// # Arguments
    ///
    /// * `mjd`, evaluation instant.
    ///
    /// # Returns
    ///
    /// `(x, y, z)` in `Kilometers` in the EclipticMeanJ2000 frame.
    #[inline]
    pub fn get_position_km(&self, mjd: ModifiedJulianDate) -> (Kilometers, Kilometers, Kilometers) {
        let offset = mjd.raw() - self.mjd_start.raw();

        if offset < Days::zero() {
            return direct_position(mjd);
        }

        let seg_idx = (offset / SEGMENT_DAYS) as usize;

        if seg_idx >= self.num_segments {
            return direct_position(mjd);
        }

        // Map jd into [-1, 1] within the segment
        let seg_start = crate::time::ModifiedJulianDate::new(
            (self.mjd_start.raw() + seg_idx as f64 * SEGMENT_DAYS).value(),
        );
        let seg_mid =
            crate::time::ModifiedJulianDate::new((seg_start.raw() + SEGMENT_DAYS * 0.5).value());
        let x = (mjd.raw() - seg_mid.raw()) / (SEGMENT_DAYS * 0.5);
        let px = Kilometers::new(cheby::evaluate(&self.cx[seg_idx], x));
        let py = Kilometers::new(cheby::evaluate(&self.cy[seg_idx], x));
        let pz = Kilometers::new(cheby::evaluate(&self.cz[seg_idx], x));
        (px, py, pz)
    }
}

#[inline]
fn direct_position(mjd: ModifiedJulianDate) -> (Kilometers, Kilometers, Kilometers) {
    let pos = DefaultEphemeris::moon_geocentric(mjd.to::<crate::JD>());
    (pos.x(), pos.y(), pos.z())
}

// =============================================================================
// NutationCache
// =============================================================================

/// Linear interpolation cache for nutation parameters (Δψ, Δε, ε₀),
/// stored as radians at regular 2-hour intervals.
pub struct NutationCache {
    /// Modified Julian Date of the first entry.
    mjd_start: ModifiedJulianDate,
    /// Number of entries.
    num_entries: usize,
    /// Pre-evaluated [dpsi_rad, deps_rad, eps0_rad] at each node.
    values: Vec<[Radians; 3]>,
}

impl NutationCache {
    /// Build the nutation cache covering `[jd_start, jd_end]`.
    ///
    /// # Arguments
    ///
    /// * `mjd_start`, inclusive lower bound of the cached window.
    /// * `mjd_end`, inclusive upper bound of the cached window.
    ///
    /// # Returns
    ///
    /// A populated [`NutationCache`] with `(Δψ, Δε, ε₀)` triplets at
    /// 2‑hour spacing across the padded window.
    pub fn new(mjd_start: ModifiedJulianDate, mjd_end: ModifiedJulianDate) -> Self {
        let pad = Days::new(1.0); // 1-day padding
        let t0 = crate::time::ModifiedJulianDate::new((mjd_start.raw() - pad).value());
        let t1 = crate::time::ModifiedJulianDate::new((mjd_end.raw() + pad).value());
        let span = t1.raw() - t0.raw();
        let num_entries = ((span / NUT_STEP_DAYS).ceil() as usize) + 1;

        let mut values = Vec::with_capacity(num_entries);
        for i in 0..num_entries {
            let mjd =
                crate::time::ModifiedJulianDate::new((t0.raw() + i as f64 * NUT_STEP_DAYS).value());
            let nut = nutation_iau2000b(mjd.to::<crate::JD>());
            values.push([nut.dpsi, nut.deps, nut.mean_obliquity]);
        }

        Self {
            mjd_start: t0,
            num_entries,
            values,
        }
    }

    /// Interpolate (Δψ, Δε, ε₀) in radians at `jd`.
    ///
    /// # Arguments
    ///
    /// * `mjd`, evaluation instant.
    ///
    /// # Returns
    ///
    /// `(dpsi, deps, eps0)` as `Radians`. Falls back to a direct
    /// IAU 2000B evaluation if `mjd` is outside the cached window.
    #[inline]
    pub fn get_nutation_rad(&self, mjd: ModifiedJulianDate) -> (Radians, Radians, Radians) {
        let offset = mjd.raw() - self.mjd_start.raw();

        if offset < Days::zero() {
            return direct_nutation(mjd);
        }

        let frac = offset / NUT_STEP_DAYS;
        let idx = frac as usize;

        if idx + 1 >= self.num_entries {
            return direct_nutation(mjd);
        }

        let t = frac - idx as f64; // interpolation parameter [0, 1)
        let a = &self.values[idx];
        let b = &self.values[idx + 1];

        (
            a[0] + t * (b[0] - a[0]),
            a[1] + t * (b[1] - a[1]),
            a[2] + t * (b[2] - a[2]),
        )
    }

    /// Build the nutation rotation matrix from cached values at `jd`.
    ///
    /// Equivalent to [`crate::astro::nutation::nutation_rotation_iau2000b`] but uses
    /// interpolated nutation parameters instead of the full 77-term series.
    ///
    /// # Arguments
    ///
    /// * `mjd`, evaluation instant.
    ///
    /// # Returns
    ///
    /// `affn::Rotation3` carrying the mean → true equator‑and‑equinox
    /// rotation `R₁(ε₀+Δε) · R₃(Δψ) · R₁(−ε₀)`.
    #[inline]
    pub fn nutation_rotation(&self, mjd: ModifiedJulianDate) -> affn::Rotation3 {
        let (dpsi, deps, eps0) = self.get_nutation_rad(mjd);

        // R1(ε0+Δε) · R3(Δψ) · R1(−ε0)
        affn::Rotation3::rx(eps0 + deps) * affn::Rotation3::rz(dpsi) * affn::Rotation3::rx(-eps0)
    }
}

#[inline]
fn direct_nutation(mjd: ModifiedJulianDate) -> (Radians, Radians, Radians) {
    let nut = nutation_iau2000b(mjd.to::<crate::JD>());
    (nut.dpsi, nut.deps, nut.mean_obliquity)
}

// =============================================================================
// MoonAltitudeContext, combines caches + site data
// =============================================================================

/// Pre-built context for fast repeated Moon altitude queries at a fixed site.
///
/// Holds:
/// * [`MoonPositionCache`], Chebyshev interpolation of ELP2000 geocentric XYZ
/// * [`NutationCache`], linear interpolation of IAU 2000B nutation values
/// * Precomputed observer ITRF position in km
///
/// The [`altitude_rad`] method reproduces the full transform chain
/// (ecliptic → equatorial → topocentric → precession → nutation → horizontal)
/// but replaces the two most expensive steps (ELP2000 and nutation) with
/// cache lookups.
pub struct MoonAltitudeContext {
    pos_cache: MoonPositionCache,
    nut_cache: NutationCache,
    /// Observer ITRF position in km: [x, y, z].
    site_itrf_km: [Kilometers; 3],
    /// Observer geodetic latitude.
    lat: crate::qtty::Radians,
    /// Observer geodetic longitude in radians (for LST computation).
    lon_rad: crate::qtty::Radians,
}

impl MoonAltitudeContext {
    /// Build an altitude context covering the MJD period for a given site.
    ///
    /// The caches are padded by 1 day on each side to accommodate Brent probes.
    ///
    /// # Arguments
    ///
    /// * `mjd_start`, inclusive lower bound of the search window.
    /// * `mjd_end`, inclusive upper bound of the search window.
    /// * `site`, geodetic observer location.
    ///
    /// # Returns
    ///
    /// A [`MoonAltitudeContext`] with both caches populated and the
    /// observer ITRF position precomputed in km.
    pub fn new(
        mjd_start: ModifiedJulianDate,
        mjd_end: ModifiedJulianDate,
        site: Geodetic<ECEF>,
    ) -> Self {
        // Convert ModifiedJulianDate to JulianDate for internal cache usage
        let pos_cache = MoonPositionCache::new(mjd_start, mjd_end);
        let nut_cache = NutationCache::new(mjd_start, mjd_end);

        // Precompute site ITRF position in km
        let site_ecef = site.to_cartesian::<Kilometer>();
        let site_itrf_km = [site_ecef.x(), site_ecef.y(), site_ecef.z()];

        Self {
            pos_cache,
            nut_cache,
            site_itrf_km,
            lat: site.lat.to::<Radian>(),
            lon_rad: site.lon.to::<Radian>(),
        }
    }

    /// Compute the Moon's topocentric altitude in radians at `mjd`.
    ///
    /// Replicates the full transform chain of
    /// [`Moon::get_horizontal`] → [`Moon::get_apparent_topocentric_equ`]
    /// but replaces ELP2000 and nutation with cached interpolations.
    ///
    /// # Arguments
    ///
    /// * `mjd`, evaluation instant.
    ///
    /// # Returns
    ///
    /// Topocentric altitude in `Radians`, no atmospheric refraction.
    #[inline]
    pub fn altitude_rad(&self, mjd: ModifiedJulianDate) -> Quantity<Radian> {
        let jd: JulianDate = mjd.to::<crate::JD>();
        let ctx: AstroContext = AstroContext::default();

        // ---------------------------------------------------------------
        // 1. Geocentric ecliptic Cartesian (km), from Chebyshev cache
        // ---------------------------------------------------------------
        let (x_ecl, y_ecl, z_ecl) = self.pos_cache.get_position_km(mjd);

        // ---------------------------------------------------------------
        // 2. EclipticMeanJ2000 → EquatorialMeanJ2000 (constant rotation about +X)
        // ---------------------------------------------------------------
        let (sin_e, cos_e) = J2000_OBLIQUITY_RAD.sin_cos();
        let x_eq = x_ecl;
        let y_eq = cos_e * y_ecl - sin_e * z_ecl;
        let z_eq = sin_e * y_ecl + cos_e * z_ecl;

        // ---------------------------------------------------------------
        // 3. Topocentric correction: subtract observer position in J2000 eq
        // ---------------------------------------------------------------
        let sx = self.site_itrf_km[0];
        let sy = self.site_itrf_km[1];
        let sz = self.site_itrf_km[2];

        // ITRF → EquatorialMeanJ2000 via full IAU 2006 + EOP chain.
        let rot_itrs = itrs_to_equatorial_mean_j2000_rotation::<
            DefaultEphemeris,
            DefaultEop,
            DefaultNutationModel,
        >(jd, &ctx);
        let [site_eq_x, site_eq_y, site_eq_z] = rot_itrs * [sx, sy, sz];

        let x_topo = x_eq - site_eq_x;
        let y_topo = y_eq - site_eq_y;
        let z_topo = z_eq - site_eq_z;

        // ---------------------------------------------------------------
        // 4. Precession: J2000 → mean-of-date
        // ---------------------------------------------------------------
        let rot_prec = precession_matrix_iau2006(mjd.to::<crate::JD>());
        let [x_mod, y_mod, z_mod] =
            rot_prec.apply_array([x_topo.value(), y_topo.value(), z_topo.value()]);

        // ---------------------------------------------------------------
        // 5. Nutation: mean-of-date → true-of-date (from cache)
        // ---------------------------------------------------------------
        let (dpsi, _deps, eps0) = self.nut_cache.get_nutation_rad(mjd);
        let rot_nut = self.nut_cache.nutation_rotation(mjd);
        let [x_tod, y_tod, z_tod] = rot_nut.apply_array([x_mod, y_mod, z_mod]);

        // ---------------------------------------------------------------
        // 6. Equatorial true-of-date → RA, Dec
        // ---------------------------------------------------------------
        let ra_rad = y_tod.atan2(x_tod);
        let r_xy = (x_tod * x_tod + y_tod * y_tod).sqrt();
        let dec_rad = z_tod.atan2(r_xy);

        // ---------------------------------------------------------------
        // 7. GAST → LST → HA → altitude
        // ---------------------------------------------------------------
        let eop = ctx.eop_at_tt(jd);
        let jd_ut1 = jd_ut1_from_tt_eop(jd, &eop);
        let gast = gast_iau2006(jd_ut1, jd, dpsi, eps0);
        let lst_rad = gast + self.lon_rad;
        let ha = (lst_rad.value() - ra_rad).rem_euclid(std::f64::consts::TAU);

        let sin_alt = dec_rad.sin() * self.lat.sin() + dec_rad.cos() * self.lat.cos() * ha.cos();

        Quantity::<Radian>::new(sin_alt.asin())
    }
}

// =============================================================================
// find_and_label_crossings, avoids probe evaluations
// =============================================================================

use crate::event::search::intervals::LabeledCrossing;
use crate::event::search::root_finding;
use crate::time::{Interval, ModifiedJulianDate};

#[allow(dead_code)]
type Mjd = ModifiedJulianDate;
type Days = crate::qtty::Quantity<crate::qtty::Day>;

/// Tiny epsilon for deduplication (same as intervals.rs).
#[allow(dead_code)]
const DEDUPE_EPS: Days = Days::new(1e-8);

/// Combined scan + Brent root-finding + labelling in a single pass.
///
/// Unlike the generic `find_crossings` + `label_crossings` pipeline in
/// `intervals.rs`, this function records the crossing direction directly
/// from the sign change that triggered the Brent solve, **eliminating
/// the 2 extra probe evaluations per crossing**.
///
/// # Arguments
///
/// * `period`, MJD search window.
/// * `step`, scan step size used to bracket sign changes.
/// * `f`, scalar function `mjd → Quantity<V>` to root‑find against `threshold`.
/// * `threshold`, value subtracted from `f` before root‑finding.
///
/// # Returns
///
/// A tuple `(crossings, start_above)` where `crossings` is the sorted,
/// deduplicated `Vec<LabeledCrossing>` of root MJDs annotated with the
/// crossing direction (`+1` entering, `−1` leaving), and `start_above`
/// is `true` if `f(start) > threshold`.
#[allow(dead_code)]
pub fn find_and_label_crossings<V, F>(
    period: Interval<ModifiedJulianDate>,
    step: Days,
    f: &F,
    threshold: crate::qtty::Quantity<V>,
) -> (Vec<LabeledCrossing>, bool)
where
    V: crate::qtty::Unit,
    F: Fn(ModifiedJulianDate) -> crate::qtty::Quantity<V>,
{
    let g = |t: Mjd| -> crate::qtty::Quantity<V> { f(t) - threshold };

    let t_start = period.start;
    let t_end = period.end;

    let start_val = g(t_start);
    let start_above = start_val > crate::qtty::Quantity::<V>::new(0.0);

    let mut labeled = Vec::new();
    let mut t = t_start;
    let mut prev = start_val;

    while t < t_end {
        let next_t = {
            let t_nx = crate::time::ModifiedJulianDate::new((t.raw() + step).value());
            if t_nx.raw() <= t_end.raw() {
                t_nx
            } else {
                t_end
            }
        };
        let next_v = g(next_t);

        if prev.signum() * next_v.signum() < 0.0 {
            if let Some(root) =
                root_finding::brent_with_values(Interval::new(t, next_t), prev, next_v, g)
            {
                if root >= t_start && root <= t_end {
                    // Direction from sign change: prev < 0 → next > 0 means entering (+1)
                    let direction = if prev < crate::qtty::Quantity::<V>::new(0.0) {
                        1
                    } else {
                        -1
                    };
                    labeled.push(LabeledCrossing { t: root, direction });
                }
            }
        }

        t = next_t;
        prev = next_v;
    }

    // Sort and deduplicate (should already be sorted from linear scan)
    labeled.sort_by(|a, b| a.t.partial_cmp(&b.t).unwrap());
    labeled.dedup_by(|a, b| {
        let dup = (a.t.raw() - b.t.raw()).abs() < DEDUPE_EPS;
        if dup {
            // Keep the earlier one (b), discard a
        }
        dup
    });

    (labeled, start_above)
}

// =============================================================================
// Tests
// =============================================================================

#[cfg(test)]
mod tests {
    use super::*;
    use crate::catalogs::observatories::ROQUE_DE_LOS_MUCHACHOS;
    use crate::event::lunar::moon_altitude_rad;
    use crate::qtty::Radians;

    #[test]
    fn chebyshev_position_accuracy() {
        // Compare cached vs. direct ELP2000 at random times within a segment
        let mjd_start: ModifiedJulianDate = crate::J2000.to::<crate::MJD>();
        let mjd_end =
            crate::time::ModifiedJulianDate::new((mjd_start.raw() + Days::new(30.0)).value());
        let cache = MoonPositionCache::new(mjd_start, mjd_end);

        for i in 0..100 {
            let mjd = crate::time::ModifiedJulianDate::new(
                (mjd_start.raw() + Days::new((i as f64) * 0.3 + 0.1)).value(),
            ); // sample every ~7 hours
            let (cx, cy, cz) = cache.get_position_km(mjd);
            let direct = DefaultEphemeris::moon_geocentric(mjd.to::<crate::JD>());
            let (dx, dy, dz) = (direct.x(), direct.y(), direct.z());
            let err = (cx - dx).abs().max((cy - dy).abs()).max((cz - dz).abs());
            // Error should be < 1 km (≈ 0.5 arcsecond at Moon distance)
            assert!(
                err < Kilometers::new(1.0),
                "Chebyshev max-axis error at MJD {mjd}: {err} (x:{} vs {}, y:{} vs {}, z:{} vs {})",
                cx,
                dx,
                cy,
                dy,
                cz,
                dz
            );
        }
    }

    #[test]
    fn moon_position_cache_falls_back_before_start() {
        let start = ModifiedJulianDate::new(60000.0);
        let end = ModifiedJulianDate::new(60010.0);
        let cache = MoonPositionCache::new(start, end);

        let t = ModifiedJulianDate::new(59998.0);
        let cached = cache.get_position_km(t);
        let direct = DefaultEphemeris::moon_geocentric(t.to::<crate::JD>());

        assert!((cached.0.value() - direct.x().value()).abs() < 1e-12);
        assert!((cached.1.value() - direct.y().value()).abs() < 1e-12);
        assert!((cached.2.value() - direct.z().value()).abs() < 1e-12);
    }

    #[test]
    fn nutation_cache_accuracy() {
        let mjd_start: ModifiedJulianDate = crate::J2000.to::<crate::MJD>();
        let mjd_end: ModifiedJulianDate =
            crate::time::ModifiedJulianDate::new((mjd_start.raw() + Days::new(30.0)).value());
        let cache = NutationCache::new(mjd_start, mjd_end);

        for i in 0..100 {
            let mjd = crate::time::ModifiedJulianDate::new(
                (mjd_start.raw() + Days::new((i as f64) * 0.3 + 0.05)).value(),
            );
            let (dpsi, deps, eps0) = cache.get_nutation_rad(mjd);
            let direct = nutation_iau2000b(mjd.to::<crate::JD>());
            let d_dpsi = direct.dpsi;
            let d_deps = direct.deps;
            let d_eps0 = direct.mean_obliquity;

            let err_dpsi = (dpsi - d_dpsi).abs();
            let err_deps = (deps - d_deps).abs();
            let err_eps0 = (eps0 - d_eps0).abs();

            // Errors should be < 5e-10 radians (≈ 0.1 milliarcseconds)
            assert!(
                err_dpsi < Radians::new(5e-10),
                "Nutation dpsi error at MJD {mjd}: {err_dpsi}"
            );
            assert!(
                err_deps < Radians::new(5e-10),
                "Nutation deps error at MJD {mjd}: {err_deps}"
            );
            assert!(
                err_eps0 < Radians::new(5e-10),
                "Nutation eps0 error at MJD {mjd}: {err_eps0}"
            );
        }
    }

    #[test]
    fn nutation_cache_falls_back_before_start() {
        let start = ModifiedJulianDate::new(60000.0);
        let end = ModifiedJulianDate::new(60010.0);
        let cache = NutationCache::new(start, end);

        let t = ModifiedJulianDate::new(59998.0);
        let cached = cache.get_nutation_rad(t);
        let direct = nutation_iau2000b(t.to::<crate::JD>());

        assert!((cached.0.value() - direct.dpsi.value()).abs() < 1e-15);
        assert!((cached.1.value() - direct.deps.value()).abs() < 1e-15);
        assert!((cached.2.value() - direct.mean_obliquity.value()).abs() < 1e-15);
    }

    #[test]
    fn cached_altitude_matches_direct() {
        let site = ROQUE_DE_LOS_MUCHACHOS.geodetic();
        let mjd_start: ModifiedJulianDate = crate::J2000.to::<crate::MJD>();
        let mjd_end =
            crate::time::ModifiedJulianDate::new((mjd_start.raw() + Days::new(7.0)).value());
        let ctx = MoonAltitudeContext::new(mjd_start, mjd_end, site);

        for i in 0..50 {
            let mjd = crate::time::ModifiedJulianDate::new(
                (mjd_start.raw() + Days::new((i as f64) * 0.14 + 0.01)).value(),
            );
            let cached_alt = ctx.altitude_rad(mjd);
            let direct_alt = moon_altitude_rad(mjd, &site);

            let err_deg = (cached_alt - direct_alt).abs().to::<Degree>();
            // Should match within ~0.01° (limited by interpolation + nutation cache)
            assert!(
                err_deg < Degrees::new(0.02),
                "Altitude error at MJD {}: cached={} direct={} err={}",
                mjd,
                cached_alt.to::<Degree>(),
                direct_alt.to::<Degree>(),
                err_deg
            );
        }
    }

    #[test]
    fn find_and_label_crossings_sine_wave() {
        // Test with a known sine wave: sin(2π(t+0.05)) crosses 0 at known times
        let f =
            |t: Mjd| Radians::new((2.0 * std::f64::consts::PI * (t.raw().value() + 0.05)).sin());
        let period = Interval::new(
            crate::time::ModifiedJulianDate::new((Days::new(0.0)).value()),
            crate::time::ModifiedJulianDate::new((Days::new(1.0)).value()),
        );
        let step = Days::new(0.01);

        let (labeled, _start_above) = find_and_label_crossings(period, step, &f, Radians::new(0.0));

        // Should find 2 crossings (at t ≈ -0.05 + 0.5 = 0.45 and t ≈ -0.05 + 1.0 = 0.95)
        assert_eq!(labeled.len(), 2, "Expected 2 crossings, got {:?}", labeled);

        // First crossing should be exiting (sin going negative), second entering
        assert_eq!(labeled[0].direction, -1);
        assert_eq!(labeled[1].direction, 1);
    }
}
