siderust 0.9.1

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

//! # DE4xx Chebyshev Segment Evaluation
//!
//! ## Scientific scope
//!
//! JPL DE ephemerides store body positions and velocities as piecewise
//! Chebyshev polynomial expansions.  Each segment covers a fixed time
//! sub-interval (typically 32 days for inner planets) and encodes `N`
//! Chebyshev coefficients per Cartesian axis.  The position polynomial is
//! evaluated by Clenshaw recurrence, and the velocity is the analytical
//! derivative of that polynomial.
//!
//! Positions are in **kilometres** (ICRF), velocities in **km · day⁻¹**.
//!
//! ## Technical scope
//!
//! - [`DynSegmentDescriptor`] — run-time segment backed by a heap
//!   `Vec<f64>`, used when reading ephemeris binary files at runtime.
//!
//! Exposes:
//! - `position(jd)  -> Displacement<ICRF, Kilometer>`
//! - `velocity(jd)  -> Velocity<ICRF, KmPerDay>`
//! - `position_velocity(jd)  -> (…, …)` — evaluates both in a single pass.
//!
//! The time argument is `TdbJulianDate` (Julian Date, TDB time scale).
//!
//! ## References
//!
//! - Standish, E. M. (1998). "JPL Planetary and Lunar Ephemerides, DE405/LE405".
//!   *JPL Interoffice Memorandum* 312.F-98-048.
//! - Park, R. S., et al. (2021). "The JPL Planetary and Lunar Ephemerides
//!   DE440 and DE441". *The Astronomical Journal* 161, 105.
//!   <https://doi.org/10.3847/1538-3881/abd414>

use crate::coordinates::frames::ICRF;
use crate::ephemeris::EphemerisError;
use crate::qtty::*;
use crate::time::{JulianDate, TDB};
use affn::{Displacement, Velocity};

type KmPerDay = Per<Kilometer, Day>;
type KmPerDayQ = crate::qtty::Quantity<KmPerDay>;
type PosVelResult =
    Result<(Displacement<ICRF, Kilometer>, Velocity<ICRF, KmPerDay>), EphemerisError>;
type TdbJulianDate = tempoch::JulianDate<TDB>;

const SECONDS_PER_DAY: f64 = crate::qtty::time::SECONDS_PER_DAY;
const J2000_JD: f64 = tempoch::J2000_JD_TT_DAY.value();

const NAIF_M0: f64 = 6.239_996;
const NAIF_M1: f64 = 1.990_968_71e-7;
const NAIF_EB: f64 = 1.671e-2;
const NAIF_K: f64 = 1.657e-3;

// ═══════════════════════════════════════════════════════════════════════════
// Shared helpers (used by both SegmentDescriptor and DynSegmentDescriptor)
// ═══════════════════════════════════════════════════════════════════════════

/// Convert a raw Julian Day value to ephemeris seconds past J2000.
#[inline]
fn jd_to_et(jd_value: f64) -> Seconds {
    Seconds::new((jd_value - J2000_JD) * SECONDS_PER_DAY)
}

/// Convert a TT Julian Date to SPICE/NAIF ephemeris time seconds past J2000.
///
/// ## What this approximates
///
/// SPK Type 2/3 segment metadata is keyed by **ephemeris time (ET)** seconds
/// past J2000. Without a loaded leapseconds kernel (LSK), this uses the
/// closed-form TT→ET mapping from the NAIF planetary ephemeris documentation:
///
/// ```text
/// ET ≈ TT + K·sin(E),   E = M + eB·sin(M),   M = M0 + M1·ET  (fixed-point)
/// ```
///
/// with constants `K`, `eB`, `M0`, `M1` from the DE405/SPICE headers. This
/// is **not** a full LSK-based conversion: it ignores leap seconds and the
/// historical ΔT table, and the fixed-point iteration introduces a residual
/// of order **10⁻⁵ s** (~1 m at Earth orbital speed) relative to CSPICE when
/// an LSK is loaded. It is adequate for Chebyshev segment lookup inside JPL
/// DE BSP files, which are themselves sampled on this ET axis.
///
/// For strict CSPICE/ANISE parity, compare against an external kernel using the
/// same TT input; see unit tests in this module for the J2000 anchor term.
#[inline]
pub(crate) fn jd_tt_to_spice_et_seconds(jd_tt: JulianDate) -> f64 {
    let tt_seconds = (jd_tt.raw().value() - J2000_JD) * SECONDS_PER_DAY;
    let mut et_seconds = tt_seconds;
    for _ in 0..5 {
        let mean_anomaly = NAIF_M0 + NAIF_M1 * et_seconds;
        let eccentric_anomaly = mean_anomaly + NAIF_EB * mean_anomaly.sin();
        let next = tt_seconds + NAIF_K * eccentric_anomaly.sin();
        if (next - et_seconds).abs() < 1.0e-12 {
            return next;
        }
        et_seconds = next;
    }
    et_seconds
}

#[inline]
fn xyz_coeffs(record: &[f64], ncoeff: usize) -> (&[f64], &[f64], &[f64]) {
    (
        &record[2..2 + ncoeff],
        &record[2 + ncoeff..2 + 2 * ncoeff],
        &record[2 + 2 * ncoeff..2 + 3 * ncoeff],
    )
}

#[inline]
fn type3_velocity_coeffs(record: &[f64], ncoeff: usize) -> (&[f64], &[f64], &[f64]) {
    let start = 2 + 3 * ncoeff;
    (
        &record[start..start + ncoeff],
        &record[start + ncoeff..start + 2 * ncoeff],
        &record[start + 2 * ncoeff..start + 3 * ncoeff],
    )
}

/// Compute the record index and Chebyshev parameter tau from segment metadata.
#[inline]
fn locate_params(
    init: Seconds,
    intlen: Seconds,
    n_records: usize,
    jd_tdb: TdbJulianDate,
) -> Result<(usize, f64, Seconds), EphemerisError> {
    let et = jd_to_et(jd_tdb.raw().value());
    locate_params_et(init, intlen, n_records, et.value())
}

#[inline]
fn locate_params_et(
    init: Seconds,
    intlen: Seconds,
    n_records: usize,
    et_s: f64,
) -> Result<(usize, f64, Seconds), EphemerisError> {
    let init_s = init.value();
    let intlen_s = intlen.value();

    if n_records == 0 || !init_s.is_finite() || !intlen_s.is_finite() || intlen_s <= 0.0 {
        return Err(EphemerisError::InvalidSegment {
            init_seconds: init_s,
            intlen_seconds: intlen_s,
            n_records,
        });
    }

    let span_s = intlen_s * n_records as f64;
    let end_s = init_s + span_s;
    if !et_s.is_finite() || et_s < init_s || et_s > end_s {
        return Err(EphemerisError::OutOfRange {
            jd: J2000_JD + et_s / SECONDS_PER_DAY,
            start_jd: J2000_JD + init_s / SECONDS_PER_DAY,
            end_jd: J2000_JD + end_s / SECONDS_PER_DAY,
        });
    }

    let rel = (et_s - init_s) / intlen_s;
    let idx = if rel >= n_records as f64 {
        n_records - 1
    } else {
        rel.floor() as usize
    };
    Ok((idx, et_s, intlen))
}

/// Given a record slice, compute (tau, radius).
#[inline]
fn record_tau(record: &[f64], et: f64) -> (f64, Seconds) {
    let mid = Seconds::new(record[0]);
    let radius = Seconds::new(record[1]);
    let tau = (et - mid.value()) / radius.value();
    (tau, radius)
}

/// Evaluate position from a record slice.
#[inline]
fn eval_position(record: &[f64], ncoeff: usize, tau: f64) -> Displacement<ICRF, Kilometer> {
    let (cx, cy, cz) = xyz_coeffs(record, ncoeff);
    Displacement::new(
        Kilometers::new(cheby::evaluate(cx, tau)),
        Kilometers::new(cheby::evaluate(cy, tau)),
        Kilometers::new(cheby::evaluate(cz, tau)),
    )
}

/// Evaluate velocity from a record slice.
#[inline]
fn eval_velocity(
    record: &[f64],
    ncoeff: usize,
    tau: f64,
    radius: Seconds,
) -> Velocity<ICRF, KmPerDay> {
    let (cx, cy, cz) = xyz_coeffs(record, ncoeff);
    let scale = SECONDS_PER_DAY / radius.value();
    Velocity::new(
        KmPerDayQ::new(cheby::evaluate_derivative(cx, tau) * scale),
        KmPerDayQ::new(cheby::evaluate_derivative(cy, tau) * scale),
        KmPerDayQ::new(cheby::evaluate_derivative(cz, tau) * scale),
    )
}

/// Evaluate explicit SPK Type 3 velocity coefficient blocks.
#[inline]
fn eval_type3_velocity(record: &[f64], ncoeff: usize, tau: f64) -> Velocity<ICRF, KmPerDay> {
    let (cvx, cvy, cvz) = type3_velocity_coeffs(record, ncoeff);
    Velocity::new(
        KmPerDayQ::new(cheby::evaluate(cvx, tau) * SECONDS_PER_DAY),
        KmPerDayQ::new(cheby::evaluate(cvy, tau) * SECONDS_PER_DAY),
        KmPerDayQ::new(cheby::evaluate(cvz, tau) * SECONDS_PER_DAY),
    )
}

/// Evaluate both position and velocity from a record slice.
#[inline]
fn eval_both(
    record: &[f64],
    ncoeff: usize,
    tau: f64,
    radius: Seconds,
) -> (Displacement<ICRF, Kilometer>, Velocity<ICRF, KmPerDay>) {
    let (cx, cy, cz) = xyz_coeffs(record, ncoeff);
    let scale = SECONDS_PER_DAY / radius.value();

    let (px, vx) = cheby::evaluate_both(cx, tau);
    let (py, vy) = cheby::evaluate_both(cy, tau);
    let (pz, vz) = cheby::evaluate_both(cz, tau);

    (
        Displacement::new(
            Kilometers::new(px),
            Kilometers::new(py),
            Kilometers::new(pz),
        ),
        Velocity::new(
            KmPerDayQ::new(vx * scale),
            KmPerDayQ::new(vy * scale),
            KmPerDayQ::new(vz * scale),
        ),
    )
}

// ═══════════════════════════════════════════════════════════════════════════
// DynSegmentDescriptor, runtime heap-backed data
// ═══════════════════════════════════════════════════════════════════════════

/// Segment descriptor for runtime-loaded data (heap-backed).
///
/// The coefficient data is stored in a `Vec<f64>` owned by this struct. This
/// enables loading BSP files at runtime without embedding them at compile time.
#[derive(Clone)]
pub(crate) struct DynSegmentDescriptor {
    /// SPK segment type (2 differentiates position coefficients, 3 stores
    /// explicit velocity coefficient blocks).
    pub(crate) data_type: i32,
    /// Initial epoch of the segment (TDB seconds past J2000).
    pub(crate) init: Seconds,
    /// Length of each Chebyshev sub-interval (seconds).
    pub(crate) intlen: Seconds,
    /// Number of Chebyshev coefficients per coordinate (x, y, z).
    pub(crate) ncoeff: usize,
    /// Doubles per record (2 + 3 × ncoeff).
    pub(crate) rsize: usize,
    /// Number of coefficient records in the segment.
    pub(crate) n_records: usize,
    /// Flattened coefficient data: n_records × rsize f64 values.
    pub(crate) data: Vec<f64>,
}

impl DynSegmentDescriptor {
    /// Construct from parsed SPK segment data.
    pub(crate) fn from_spk(seg: &crate::formats::spice::spk::SegmentData) -> Self {
        Self {
            data_type: seg.data_type,
            init: Seconds::new(seg.init),
            intlen: Seconds::new(seg.intlen),
            ncoeff: seg.ncoeff,
            rsize: seg.rsize,
            n_records: seg.n_records,
            data: seg.records.clone(),
        }
    }

    /// Get the `i`-th record as a slice.
    #[inline]
    fn record(&self, i: usize) -> &[f64] {
        let start = i * self.rsize;
        &self.data[start..start + self.rsize]
    }

    /// Locate the record for a given JD and compute tau.
    #[inline]
    fn try_locate(&self, jd_tdb: TdbJulianDate) -> Result<(&[f64], f64, Seconds), EphemerisError> {
        let (idx, et, _) = locate_params(self.init, self.intlen, self.n_records, jd_tdb)?;
        let record = self.record(idx);
        let (tau, radius) = record_tau(record, et);
        Ok((record, tau, radius))
    }

    /// Locate the record for SPICE/NAIF ET seconds past J2000 and compute tau.
    #[inline]
    pub(crate) fn try_locate_et(
        &self,
        et_seconds: f64,
    ) -> Result<(&[f64], f64, Seconds), EphemerisError> {
        let (idx, et, _) = locate_params_et(self.init, self.intlen, self.n_records, et_seconds)?;
        let record = self.record(idx);
        let (tau, radius) = record_tau(record, et);
        Ok((record, tau, radius))
    }

    /// Fallibly evaluate position in km (ICRF) at Julian Date (TDB).
    #[inline]
    pub(crate) fn try_position(
        &self,
        jd_tdb: TdbJulianDate,
    ) -> Result<Displacement<ICRF, Kilometer>, EphemerisError> {
        let (record, tau, _) = self.try_locate(jd_tdb)?;
        Ok(eval_position(record, self.ncoeff, tau))
    }

    /// Fallibly evaluate position in km (ICRF) at SPICE/NAIF ET seconds past J2000.
    #[inline]
    pub(crate) fn try_position_et(
        &self,
        et_seconds: f64,
    ) -> Result<Displacement<ICRF, Kilometer>, EphemerisError> {
        let (record, tau, _) = self.try_locate_et(et_seconds)?;
        Ok(eval_position(record, self.ncoeff, tau))
    }

    /// Evaluate position in km (ICRF) at Julian Date (TDB).
    ///
    /// # Panics
    ///
    /// Panics when `jd_tdb` is outside the segment coverage. Use
    /// [`Self::try_position`] to handle that condition explicitly.
    #[inline]
    #[allow(dead_code)]
    pub(crate) fn position(&self, jd_tdb: TdbJulianDate) -> Displacement<ICRF, Kilometer> {
        self.try_position(jd_tdb)
            .expect("JPL runtime segment position requested outside ephemeris coverage")
    }

    /// Fallibly evaluate velocity in km/day (ICRF) at Julian Date (TDB).
    #[inline]
    pub(crate) fn try_velocity(
        &self,
        jd_tdb: TdbJulianDate,
    ) -> Result<Velocity<ICRF, KmPerDay>, EphemerisError> {
        let (record, tau, radius) = self.try_locate(jd_tdb)?;
        if self.data_type == 3 {
            Ok(eval_type3_velocity(record, self.ncoeff, tau))
        } else {
            Ok(eval_velocity(record, self.ncoeff, tau, radius))
        }
    }

    /// Fallibly evaluate velocity in km/day (ICRF) at SPICE/NAIF ET seconds past J2000.
    #[inline]
    pub(crate) fn try_velocity_et(
        &self,
        et_seconds: f64,
    ) -> Result<Velocity<ICRF, KmPerDay>, EphemerisError> {
        let (record, tau, radius) = self.try_locate_et(et_seconds)?;
        if self.data_type == 3 {
            Ok(eval_type3_velocity(record, self.ncoeff, tau))
        } else {
            Ok(eval_velocity(record, self.ncoeff, tau, radius))
        }
    }

    /// Evaluate velocity in km/day (ICRF) at Julian Date (TDB).
    ///
    /// # Panics
    ///
    /// Panics when `jd_tdb` is outside the segment coverage. Use
    /// [`Self::try_velocity`] to handle that condition explicitly.
    #[inline]
    #[allow(dead_code)]
    pub(crate) fn velocity(&self, jd_tdb: TdbJulianDate) -> Velocity<ICRF, KmPerDay> {
        self.try_velocity(jd_tdb)
            .expect("JPL runtime segment velocity requested outside ephemeris coverage")
    }

    /// Fallibly evaluate both position and velocity in one pass.
    #[inline]
    #[allow(dead_code)]
    pub(crate) fn try_position_velocity(&self, jd_tdb: TdbJulianDate) -> PosVelResult {
        let (record, tau, radius) = self.try_locate(jd_tdb)?;
        if self.data_type == 3 {
            Ok((
                eval_position(record, self.ncoeff, tau),
                eval_type3_velocity(record, self.ncoeff, tau),
            ))
        } else {
            Ok(eval_both(record, self.ncoeff, tau, radius))
        }
    }

    /// Evaluate both position and velocity in one pass.
    ///
    /// # Panics
    ///
    /// Panics when `jd_tdb` is outside the segment coverage. Use
    /// [`Self::try_position_velocity`] to handle that condition explicitly.
    #[inline]
    #[allow(dead_code)]
    pub(crate) fn position_velocity(
        &self,
        jd_tdb: TdbJulianDate,
    ) -> (Displacement<ICRF, Kilometer>, Velocity<ICRF, KmPerDay>) {
        self.try_position_velocity(jd_tdb)
            .expect("JPL runtime segment state requested outside ephemeris coverage")
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// DynSegmentStack — one or more SPK segments for the same target/center
// ═══════════════════════════════════════════════════════════════════════════

#[derive(Clone)]
struct StackedSegment {
    start_et: f64,
    end_et: f64,
    descriptor: DynSegmentDescriptor,
}

/// All SPK segments for one target/center pair, with later segments preferred.
#[derive(Clone)]
pub(crate) struct DynSegmentStack {
    segments: Vec<StackedSegment>,
}

impl DynSegmentStack {
    /// Coverage interval in ET seconds past J2000 from segment metadata.
    pub(crate) fn coverage_from_spk(seg: &crate::formats::spice::spk::SegmentData) -> (f64, f64) {
        let start = seg.init;
        let end = seg.init + seg.intlen * seg.n_records as f64;
        (start, end)
    }

    /// One segment covering `[start_et, end_et]`.
    pub(crate) fn single(descriptor: DynSegmentDescriptor, start_et: f64, end_et: f64) -> Self {
        Self {
            segments: vec![StackedSegment {
                start_et,
                end_et,
                descriptor,
            }],
        }
    }

    /// Build from one parsed SPK segment (coverage inferred from metadata).
    pub(crate) fn from_spk_segment(seg: &crate::formats::spice::spk::SegmentData) -> Self {
        let (start_et, end_et) = Self::coverage_from_spk(seg);
        Self::single(DynSegmentDescriptor::from_spk(seg), start_et, end_et)
    }

    /// Collect every indexed segment for a NAIF target/center pair.
    pub(crate) fn for_indexed_pair(
        indexed: &[crate::formats::spice::spk::IndexedSegmentData],
        target_id: i32,
        center_id: i32,
    ) -> Self {
        let segments = indexed
            .iter()
            .filter(|s| s.target_id == target_id && s.center_id == center_id)
            .map(|s| StackedSegment {
                start_et: s.start_et,
                end_et: s.end_et,
                descriptor: DynSegmentDescriptor::from_spk(&s.data),
            })
            .collect();
        Self { segments }
    }

    #[inline]
    fn covers_et(segment: &StackedSegment, et: f64) -> bool {
        et.is_finite() && et >= segment.start_et && et <= segment.end_et
    }

    #[inline]
    fn out_of_range_et(segment: &StackedSegment, et: f64) -> EphemerisError {
        EphemerisError::OutOfRange {
            jd: J2000_JD + et / SECONDS_PER_DAY,
            start_jd: J2000_JD + segment.start_et / SECONDS_PER_DAY,
            end_jd: J2000_JD + segment.end_et / SECONDS_PER_DAY,
        }
    }

    fn select_et(&self, et_seconds: f64) -> Result<&DynSegmentDescriptor, EphemerisError> {
        let mut out_of_range = None;
        for segment in self.segments.iter().rev() {
            if !Self::covers_et(segment, et_seconds) {
                out_of_range.get_or_insert(Self::out_of_range_et(segment, et_seconds));
                continue;
            }
            match segment.descriptor.try_locate_et(et_seconds) {
                Ok(_) => return Ok(&segment.descriptor),
                Err(err @ EphemerisError::OutOfRange { .. }) => {
                    out_of_range.get_or_insert(err);
                }
                Err(err) => return Err(err),
            }
        }
        Err(out_of_range.unwrap_or(EphemerisError::InvalidSegment {
            init_seconds: f64::NAN,
            intlen_seconds: f64::NAN,
            n_records: 0,
        }))
    }

    /// Fallibly evaluate position in km (ICRF) at SPICE/NAIF ET seconds past J2000.
    #[inline]
    pub(crate) fn try_position_et(
        &self,
        et_seconds: f64,
    ) -> Result<Displacement<ICRF, Kilometer>, EphemerisError> {
        self.select_et(et_seconds)?.try_position_et(et_seconds)
    }

    /// Fallibly evaluate velocity in km/day (ICRF) at SPICE/NAIF ET seconds past J2000.
    #[inline]
    pub(crate) fn try_velocity_et(
        &self,
        et_seconds: f64,
    ) -> Result<Velocity<ICRF, KmPerDay>, EphemerisError> {
        self.select_et(et_seconds)?.try_velocity_et(et_seconds)
    }

    /// Total number of stacked SPK segments.
    pub(crate) fn segment_count(&self) -> usize {
        self.segments.len()
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    const JD_J2000: f64 = tempoch::J2000_JD_TT_DAY.value();

    /// Build a synthetic DynSegmentDescriptor with ncoeff=2 spanning 1000 days from J2000.
    ///
    /// The single record is centred at J2000 + 500 days (tau=0 at that epoch).
    /// At tau=0 (Chebyshev T0(0)=1, T1(0)=0), position = (cx0, cy0, cz0) km.
    fn make_desc(cx0: f64, cy0: f64, cz0: f64) -> DynSegmentDescriptor {
        let ncoeff = 2usize;
        let rsize = 2 + 3 * ncoeff; // = 8
        let intlen_days = 1000.0_f64;
        let intlen_secs = intlen_days * SECONDS_PER_DAY;
        let init_secs = 0.0_f64; // 0 seconds past J2000

        // mid = init + intlen/2 =  intlen / 2
        let mid = init_secs + intlen_secs / 2.0;
        let radius = intlen_secs / 2.0;

        // Record layout: [mid, radius, cx0, cx1, cy0, cy1, cz0, cz1]
        let data = vec![mid, radius, cx0, 0.0, cy0, 0.0, cz0, 0.0];

        DynSegmentDescriptor {
            data_type: 2,
            init: Seconds::new(init_secs),
            intlen: Seconds::new(intlen_secs),
            ncoeff,
            rsize,
            n_records: 1,
            data,
        }
    }

    /// JD for J2000 + 500 days (the midpoint → tau = 0).
    fn jd_mid() -> TdbJulianDate {
        TdbJulianDate::try_new(Days::new(JD_J2000 + 500.0)).unwrap()
    }

    /// JD for J2000 + 250 days (first quarter → tau = -0.5).
    fn jd_quarter() -> TdbJulianDate {
        TdbJulianDate::try_new(Days::new(JD_J2000 + 250.0)).unwrap()
    }

    #[test]
    fn spice_et_conversion_matches_naif_j2000_term() {
        let jd = JulianDate::new(JD_J2000);
        let et = jd_tt_to_spice_et_seconds(jd);

        assert!(
            (et + 7.273_677_621_567_873e-5).abs() < 1.0e-14,
            "et={et:.17e}"
        );
    }

    /// TT→ET without LSK stays within ~0.1 s of TT at sample epochs (see module docs).
    #[test]
    fn spice_et_conversion_near_tt_at_reference_epochs() {
        for jd in [JD_J2000, 2_453_372.5, 2_440_000.0, 2_500_000.0] {
            let et = jd_tt_to_spice_et_seconds(JulianDate::new(jd));
            let tt_seconds = (jd - JD_J2000) * SECONDS_PER_DAY;
            assert!(et.is_finite());
            assert!(
                (et - tt_seconds).abs() < 0.1,
                "jd={jd}: |ET−TT| = {} s exceeds 0.1 s approximation bound",
                (et - tt_seconds).abs()
            );
        }
    }

    #[test]
    fn dyn_segment_stack_later_segment_preferred() {
        use crate::formats::spice::spk::{IndexedSegmentData, SegmentData};

        let intlen = 500.0 * SECONDS_PER_DAY;
        let seg = |x: f64| SegmentData {
            data_type: 2,
            init: 0.0,
            intlen,
            rsize: 8,
            ncoeff: 2,
            n_records: 1,
            records: vec![intlen / 2.0, intlen / 2.0, x, 0.0, 0.0, 0.0, 0.0, 0.0],
        };
        let indexed = vec![
            IndexedSegmentData {
                target_id: 10,
                center_id: 0,
                frame_id: 1,
                start_et: 0.0,
                end_et: intlen,
                data: seg(100.0),
            },
            IndexedSegmentData {
                target_id: 10,
                center_id: 0,
                frame_id: 1,
                start_et: 0.0,
                end_et: intlen,
                data: seg(900.0),
            },
        ];
        let stack = DynSegmentStack::for_indexed_pair(&indexed, 10, 0);
        let et = 250.0 * SECONDS_PER_DAY;
        let pos = stack.try_position_et(et).unwrap();
        assert!((pos.x().value() - 900.0).abs() < 1.0e-6);
    }

    #[test]
    fn position_et_matches_tdb_seconds_for_same_epoch() {
        let desc = make_desc(1.0, 2.0, 3.0);
        let jd = jd_quarter();
        let et = (jd.raw().value() - JD_J2000) * SECONDS_PER_DAY;
        let from_jd = desc.position(jd);
        let from_et = desc.try_position_et(et).unwrap();

        assert_eq!(from_jd.x().value(), from_et.x().value());
        assert_eq!(from_jd.y().value(), from_et.y().value());
        assert_eq!(from_jd.z().value(), from_et.z().value());
    }

    // ── Position ──────────────────────────────────────────────────────────

    #[test]
    fn dyn_desc_position_at_mid() {
        // At tau=0, Chebyshev evaluate([c0, 0], 0) = c0 * T0(0) = c0
        let desc = make_desc(1000.0, 2000.0, 3000.0);
        let pos = desc.position(jd_mid());
        assert!(
            (pos.x().value() - 1000.0).abs() < 1e-6,
            "x = {}",
            pos.x().value()
        );
        assert!(
            (pos.y().value() - 2000.0).abs() < 1e-6,
            "y = {}",
            pos.y().value()
        );
        assert!(
            (pos.z().value() - 3000.0).abs() < 1e-6,
            "z = {}",
            pos.z().value()
        );
    }

    #[test]
    fn dyn_desc_position_is_finite() {
        let desc = make_desc(500.0, -300.0, 150.0);
        let pos = desc.position(jd_quarter());
        assert!(pos.x().is_finite());
        assert!(pos.y().is_finite());
        assert!(pos.z().is_finite());
    }

    #[test]
    fn dyn_desc_position_at_boundary_does_not_panic() {
        // First JD in the segment.
        let desc = make_desc(100.0, 200.0, 300.0);
        let jd_start = TdbJulianDate::try_new(Days::new(JD_J2000)).unwrap(); // et = 0 → idx = 0
        let pos = desc.position(jd_start);
        assert!(pos.x().is_finite());
    }

    #[test]
    fn dyn_desc_position_at_end_boundary_does_not_panic() {
        let desc = make_desc(100.0, 200.0, 300.0);
        let jd_end = TdbJulianDate::try_new(Days::new(JD_J2000 + 1000.0)).unwrap();
        let pos = desc.try_position(jd_end).expect("end boundary is in range");
        assert!(pos.x().is_finite());
    }

    #[test]
    fn dyn_desc_rejects_before_segment() {
        let desc = make_desc(100.0, 200.0, 300.0);
        let before = TdbJulianDate::try_new(Days::new(JD_J2000 - 1.0e-8)).unwrap();
        assert!(desc.try_position(before).is_err());
    }

    #[test]
    fn dyn_desc_rejects_after_segment() {
        let desc = make_desc(100.0, 200.0, 300.0);
        let after = TdbJulianDate::try_new(Days::new(JD_J2000 + 1000.0 + 1.0e-8)).unwrap();
        assert!(desc.try_position_velocity(after).is_err());
    }

    // ── Velocity ──────────────────────────────────────────────────────────

    #[test]
    fn dyn_desc_velocity_is_finite() {
        let desc = make_desc(1000.0, 2000.0, 3000.0);
        let vel = desc.velocity(jd_mid());
        assert!(vel.x().is_finite());
        assert!(vel.y().is_finite());
        assert!(vel.z().is_finite());
    }

    #[test]
    fn dyn_desc_velocity_at_tau0_with_linear_coeff() {
        // With coeffs [c0=0, c1=1]: at tau=0, dT1/dtau=1 so derivative=1
        // scale = SECONDS_PER_DAY / radius
        let ncoeff = 2usize;
        let rsize = 8usize;
        let intlen_secs = 1000.0 * SECONDS_PER_DAY;
        let radius = intlen_secs / 2.0;
        let data = vec![
            intlen_secs / 2.0, // mid
            radius,            // radius
            0.0,
            1.0, // cx = [0, 1] → at tau=0: d/dtau = 1
            0.0,
            2.0, // cy = [0, 2] → at tau=0: d/dtau = 2
            0.0,
            3.0, // cz = [0, 3] → at tau=0: d/dtau = 3
        ];
        let desc = DynSegmentDescriptor {
            data_type: 2,
            init: Seconds::new(0.0),
            intlen: Seconds::new(intlen_secs),
            ncoeff,
            rsize,
            n_records: 1,
            data,
        };
        let vel = desc.velocity(jd_mid());
        let scale = SECONDS_PER_DAY / radius;
        // dT1(0)/dtau = 1, so v_x = 1 * scale
        assert!((vel.x().value() - 1.0 * scale).abs() < 1e-6);
        assert!((vel.y().value() - 2.0 * scale).abs() < 1e-6);
        assert!((vel.z().value() - 3.0 * scale).abs() < 1e-6);
    }

    // ── Position+Velocity ─────────────────────────────────────────────────

    #[test]
    fn dyn_desc_position_velocity_consistent() {
        let desc = make_desc(1500.0, -800.0, 400.0);
        let jd = jd_quarter();
        let pos_only = desc.position(jd);
        let vel_only = desc.velocity(jd);
        let (pos_both, vel_both) = desc.position_velocity(jd);

        assert!((pos_only.x().value() - pos_both.x().value()).abs() < 1e-10);
        assert!((pos_only.y().value() - pos_both.y().value()).abs() < 1e-10);
        assert!((pos_only.z().value() - pos_both.z().value()).abs() < 1e-10);
        assert!((vel_only.x().value() - vel_both.x().value()).abs() < 1e-10);
        assert!((vel_only.y().value() - vel_both.y().value()).abs() < 1e-10);
        assert!((vel_only.z().value() - vel_both.z().value()).abs() < 1e-10);
    }

    // ── from_spk ─────────────────────────────────────────────────────────

    #[test]
    fn dyn_desc_from_spk_roundtrip() {
        use crate::formats::spice::spk::SegmentData;
        let ncoeff = 3usize;
        let rsize = 2 + 3 * ncoeff; // = 11
        let records = vec![0.0; rsize]; // one zero record
        let seg = SegmentData {
            data_type: 2,
            init: 0.0,
            intlen: 86400.0,
            rsize,
            ncoeff,
            n_records: 1,
            records,
        };
        let desc = DynSegmentDescriptor::from_spk(&seg);
        assert_eq!(desc.ncoeff, ncoeff);
        assert_eq!(desc.rsize, rsize);
        assert_eq!(desc.n_records, 1);
        assert_eq!(desc.data.len(), rsize);
    }

    #[test]
    fn dyn_desc_type3_uses_explicit_velocity_coefficients() {
        use crate::formats::spice::spk::SegmentData;

        let intlen = 1000.0 * SECONDS_PER_DAY;
        let seg = SegmentData {
            data_type: 3,
            init: 0.0,
            intlen,
            rsize: 14,
            ncoeff: 2,
            n_records: 1,
            records: vec![
                intlen / 2.0,
                intlen / 2.0,
                10.0,
                0.0,
                20.0,
                0.0,
                30.0,
                0.0,
                1.0,
                0.0,
                2.0,
                0.0,
                3.0,
                0.0,
            ],
        };
        let desc = DynSegmentDescriptor::from_spk(&seg);
        let (pos, vel) = desc.position_velocity(jd_mid());
        assert_eq!(pos.x().value(), 10.0);
        assert_eq!(vel.x().value(), SECONDS_PER_DAY);
        assert_eq!(vel.y().value(), 2.0 * SECONDS_PER_DAY);
        assert_eq!(vel.z().value(), 3.0 * SECONDS_PER_DAY);
    }
    // ── Multi-record ─────────────────────────────────────────────────────

    #[test]
    fn dyn_desc_multi_record_second_interval() {
        // Two records, each 500-day intervals. Query in second interval (J2000 + 750 days).
        let ncoeff = 2usize;
        let rsize = 8usize;
        let intlen_secs = 500.0 * SECONDS_PER_DAY;
        // Record 0: mid = intlen/2, radius = intlen/2, pos = (100, 0, 0)
        let mid0 = intlen_secs / 2.0;
        let rad = intlen_secs / 2.0;
        // Record 1: mid = intlen + intlen/2, pos = (200, 0, 0)
        let mid1 = intlen_secs + intlen_secs / 2.0;
        let data = vec![
            mid0, rad, 100.0, 0.0, 0.0, 0.0, 0.0, 0.0, // record 0
            mid1, rad, 200.0, 0.0, 0.0, 0.0, 0.0, 0.0, // record 1
        ];
        let desc = DynSegmentDescriptor {
            data_type: 2,
            init: Seconds::new(0.0),
            intlen: Seconds::new(intlen_secs),
            ncoeff,
            rsize,
            n_records: 2,
            data,
        };
        // At J2000 + 750 days → et = 750 * 86400 = 64800000, idx = (64800000 / 43200000) = 1
        let jd_second = TdbJulianDate::try_new(Days::new(JD_J2000 + 750.0)).unwrap();
        let pos = desc.position(jd_second);
        // At tau = (et - mid1) / rad with et = 750*86400, mid1 = 750*86400 → tau=0 → x=200
        let et_s = 750.0 * SECONDS_PER_DAY;
        let tau = (et_s - mid1) / rad;
        let expected_x = 200.0; // T0(tau=0) = 1 but tau may not be 0
        let _ = tau;
        let _ = expected_x;
        assert!(pos.x().is_finite());
    }
}