principia 0.2.0

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

//! Adaptive propagation driver, events, and results.
//!
//! ## Scientific scope
//!
//! Couples an adaptive Runge-Kutta integrator to a Cartesian acceleration
//! model while handling time bounds, output sampling, and event detection.
//!
//! ## Technical scope
//!
//! [`PropagationConfig`] carries the time window, step-size bounds,
//! tolerances, event detectors, and output schedule. [`propagate`] is generic
//! over the caller-owned context type, time scale, center, and frame.
//!
//! ## References
//!
//! * Hairer, Nørsett, Wanner, *Solving Ordinary Differential Equations I*, §II.5.

use core::fmt;

#[cfg(any(feature = "alloc", feature = "std"))]
use alloc::boxed::Box;
#[cfg(any(feature = "alloc", feature = "std"))]
use alloc::vec::Vec;

use affn::centers::ReferenceCenter;
use affn::frames::ReferenceFrame;
use qtty::length::Kilometers;
use tempoch::ContinuousScale;

#[cfg(any(feature = "alloc", feature = "std"))]
use qtty::IntegratorTolerances;
#[cfg(any(feature = "alloc", feature = "std"))]
use qtty::Second;
#[cfg(any(feature = "alloc", feature = "std"))]
use tempoch::Time;

use crate::error::PrincipiaError;
#[cfg(any(feature = "alloc", feature = "std"))]
use crate::integrators::AdaptiveStepper;
#[cfg(any(feature = "alloc", feature = "std"))]
use crate::models::AccelerationModel;
use crate::state::DynamicsState;

/// Error family produced by the propagation driver.
#[derive(Debug)]
pub enum PropagationError {
    /// The supplied propagation configuration is inconsistent.
    InvalidConfiguration(PrincipiaError),
    /// Wrap an underlying stepper failure.
    StepControl(PrincipiaError),
    /// The step controller requested a step smaller than the configured minimum.
    StepBelowMinimum {
        /// Requested step magnitude in seconds.
        h_requested: f64,
        /// Configured minimum allowed step magnitude in seconds.
        h_min: f64,
    },
    /// The driver exceeded the configured accepted-step budget.
    MaxStepsExceeded {
        /// Configured maximum number of accepted steps.
        max_steps: usize,
    },
    /// An event detector failed while evaluating its switching function.
    EventEvaluation {
        /// Detector name.
        name: &'static str,
        /// Underlying detector error.
        source: PrincipiaError,
    },
}

impl fmt::Display for PropagationError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            Self::InvalidConfiguration(source) => {
                write!(f, "invalid propagation configuration: {source}")
            }
            Self::StepControl(source) => write!(f, "integrator step control error: {source}"),
            Self::StepBelowMinimum { h_requested, h_min } => write!(
                f,
                "requested step {h_requested:e} s falls below configured minimum {h_min:e} s"
            ),
            Self::MaxStepsExceeded { max_steps } => {
                write!(f, "propagation exceeded max_steps={max_steps}")
            }
            Self::EventEvaluation { name, source } => {
                write!(f, "event '{name}' evaluation failed: {source}")
            }
        }
    }
}

impl core::error::Error for PropagationError {
    fn source(&self) -> Option<&(dyn core::error::Error + 'static)> {
        match self {
            Self::InvalidConfiguration(source) => Some(source),
            Self::StepControl(source) => Some(source),
            Self::EventEvaluation { source, .. } => Some(source),
            _ => None,
        }
    }
}

/// Direction of a state-function crossing that constitutes an event.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub enum EventDirection {
    /// Trigger on rising zero-crossings (state function goes from negative to positive).
    Rising,
    /// Trigger on falling zero-crossings (state function goes from positive to negative).
    Falling,
    /// Trigger on either direction of crossing.
    Either,
}

/// Propagation configuration for the low-level driver.
#[cfg(any(feature = "alloc", feature = "std"))]
// serde: not derived — contains Box<dyn EventDetector<...>>.
pub struct PropagationConfig<Ctx, S, C, F>
where
    S: ContinuousScale,
    C: ReferenceCenter,
    F: ReferenceFrame,
{
    /// Start epoch of the propagation window.
    pub t_start: Time<S>,
    /// End epoch of the propagation window.
    pub t_end: Time<S>,
    /// Initial step-size guess.
    pub h0: Second,
    /// Maximum allowed step magnitude.
    pub h_max: Second,
    /// Minimum allowed step magnitude.
    pub h_min: Second,
    /// Adaptive-step tolerances.
    pub tolerances: IntegratorTolerances,
    /// Optional regular output cadence.
    pub output_every: Option<Second>,
    /// Additional explicit output epochs.
    pub output_at: Vec<Time<S>>,
    /// Registered event detectors.
    pub events: Vec<Box<dyn EventDetector<Ctx, S, C, F>>>,
    /// Maximum number of accepted steps.
    pub max_steps: usize,
}

#[cfg(any(feature = "alloc", feature = "std"))]
impl<Ctx, S, C, F> PropagationConfig<Ctx, S, C, F>
where
    S: ContinuousScale,
    C: ReferenceCenter,
    F: ReferenceFrame,
{
    /// Construct a propagation configuration with conservative defaults.
    pub fn new(t_start: Time<S>, t_end: Time<S>) -> Self {
        Self {
            t_start,
            t_end,
            h0: Second::new(30.0),
            h_max: Second::new(86_400.0),
            h_min: Second::new(1.0e-6),
            tolerances: IntegratorTolerances::uniform(1e-9, 1e-6, 1e-9),
            output_every: None,
            output_at: Vec::new(),
            events: Vec::new(),
            max_steps: 1_000_000,
        }
    }

    /// Override the initial step guess.
    pub fn with_initial_step(mut self, h0: Second) -> Self {
        self.h0 = h0;
        self
    }

    /// Override the maximum step magnitude.
    pub fn with_max_step(mut self, h_max: Second) -> Self {
        self.h_max = h_max;
        self
    }

    /// Override the minimum step magnitude.
    pub fn with_min_step(mut self, h_min: Second) -> Self {
        self.h_min = h_min;
        self
    }

    /// Override the adaptive-step tolerances.
    pub fn with_tolerances(mut self, tolerances: IntegratorTolerances) -> Self {
        self.tolerances = tolerances;
        self
    }

    /// Request outputs at a regular cadence.
    pub fn with_output_every(mut self, dt: Second) -> Self {
        self.output_every = Some(dt);
        self
    }

    /// Request outputs at explicit epochs.
    pub fn with_output_at(mut self, times: Vec<Time<S>>) -> Self {
        self.output_at = times;
        self
    }

    /// Append one event detector.
    pub fn with_event(mut self, ev: Box<dyn EventDetector<Ctx, S, C, F>>) -> Self {
        self.events.push(ev);
        self
    }

    /// Override the accepted-step budget.
    pub fn with_max_steps(mut self, max_steps: usize) -> Self {
        self.max_steps = max_steps;
        self
    }

    /// Return the signed propagation duration in seconds.
    pub fn total_duration_s(&self) -> f64 {
        (self.t_end - self.t_start).value()
    }

    /// Validates this configuration, returning an error if any parameter is inconsistent.
    pub fn validate(&self) -> Result<(), PrincipiaError> {
        if !self.total_duration_s().is_finite() {
            return Err(PrincipiaError::InvalidPropagationConfig {
                reason: "PropagationConfig: propagation span must be finite",
            });
        }
        if !self.h0.value().is_finite() || self.h0.value() == 0.0 {
            return Err(PrincipiaError::InvalidPropagationConfig {
                reason: "PropagationConfig: h0 must be finite and non-zero",
            });
        }
        if !self.h_min.value().is_finite() || self.h_min.value() <= 0.0 {
            return Err(PrincipiaError::InvalidPropagationConfig {
                reason: "PropagationConfig: h_min must be finite and positive",
            });
        }
        if !self.h_max.value().is_finite() || self.h_max.value() <= 0.0 {
            return Err(PrincipiaError::InvalidPropagationConfig {
                reason: "PropagationConfig: h_max must be finite and positive",
            });
        }
        if self.h_min.value().abs() > self.h_max.value().abs() {
            return Err(PrincipiaError::InvalidPropagationConfig {
                reason: "PropagationConfig: h_min must not exceed h_max",
            });
        }
        if !self.tolerances.rel.value().is_finite() || self.tolerances.rel.value() <= 0.0 {
            return Err(PrincipiaError::InvalidPropagationConfig {
                reason: "PropagationConfig: relative tolerance must be finite and positive",
            });
        }
        for abs_tol in self.tolerances.abs_pos {
            if !abs_tol.value().is_finite() || abs_tol.value() <= 0.0 {
                return Err(PrincipiaError::InvalidPropagationConfig {
                    reason:
                        "PropagationConfig: absolute position tolerance must be finite and positive",
                });
            }
        }
        for abs_tol in self.tolerances.abs_vel {
            if !abs_tol.value().is_finite() || abs_tol.value() <= 0.0 {
                return Err(PrincipiaError::InvalidPropagationConfig {
                    reason:
                        "PropagationConfig: absolute velocity tolerance must be finite and positive",
                });
            }
        }
        if let Some(output_every) = self.output_every {
            if !output_every.value().is_finite() || output_every.value() <= 0.0 {
                return Err(PrincipiaError::InvalidPropagationConfig {
                    reason: "PropagationConfig: output_every must be finite and positive",
                });
            }
        }
        if self.max_steps == 0 {
            return Err(PrincipiaError::InvalidPropagationConfig {
                reason: "PropagationConfig: max_steps must be at least 1",
            });
        }
        Ok(())
    }
}

/// Zero-crossing event detector evaluated on propagated states.
pub trait EventDetector<Ctx, S, C, F>: Send + Sync
where
    S: ContinuousScale,
    C: ReferenceCenter,
    F: ReferenceFrame,
{
    /// Short event identifier.
    fn name(&self) -> &'static str;

    /// Evaluate the detector switching function `g(t, x)`.
    fn evaluate(&self, state: &DynamicsState<S, C, F>, ctx: &Ctx) -> Result<f64, PrincipiaError>;

    /// Return `true` if the event terminates propagation.
    fn terminal(&self) -> bool {
        false
    }

    /// Return `true` if a sign change from `g_before` to `g_after` should be
    /// treated as an event occurrence.
    fn accepts_crossing(&self, g_before: f64, g_after: f64) -> bool {
        g_before == 0.0 || g_after == 0.0 || (g_before * g_after) < 0.0
    }
}

/// Detector for a radial threshold crossing `|r| = threshold`.
#[derive(Debug, Clone, Copy)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct RadialThresholdEvent {
    /// Trigger radius.
    pub threshold: Kilometers,
    /// Accepted crossing direction.
    pub direction: EventDirection,
    /// If `true`, stop propagation when the event fires.
    pub terminal: bool,
}

impl RadialThresholdEvent {
    /// Construct a non-terminal radial-threshold detector.
    pub fn new(threshold: Kilometers) -> Self {
        Self {
            threshold,
            direction: EventDirection::Rising,
            terminal: false,
        }
    }

    /// Configure the detector as terminal or non-terminal.
    pub fn terminal(mut self, terminal: bool) -> Self {
        self.terminal = terminal;
        self
    }

    /// Configure the accepted crossing direction.
    pub fn direction(mut self, direction: EventDirection) -> Self {
        self.direction = direction;
        self
    }
}

impl<Ctx, S, C, F> EventDetector<Ctx, S, C, F> for RadialThresholdEvent
where
    S: ContinuousScale,
    C: ReferenceCenter,
    F: ReferenceFrame,
{
    fn name(&self) -> &'static str {
        "radial_threshold"
    }

    fn evaluate(&self, state: &DynamicsState<S, C, F>, _ctx: &Ctx) -> Result<f64, PrincipiaError> {
        let x = state.position.x().value();
        let y = state.position.y().value();
        let z = state.position.z().value();
        Ok((x * x + y * y + z * z).sqrt() - self.threshold.value())
    }

    fn terminal(&self) -> bool {
        self.terminal
    }

    fn accepts_crossing(&self, g_before: f64, g_after: f64) -> bool {
        match self.direction {
            EventDirection::Rising => g_before < 0.0 && g_after >= 0.0,
            EventDirection::Falling => g_before > 0.0 && g_after <= 0.0,
            EventDirection::Either => {
                g_before == 0.0 || g_after == 0.0 || (g_before * g_after) < 0.0
            }
        }
    }
}

/// Recorded event occurrence.
#[derive(Debug, Clone)]
pub struct EventOccurrence<S, C, F>
where
    S: ContinuousScale,
    C: ReferenceCenter,
    F: ReferenceFrame,
{
    /// Event identifier.
    pub event_name: &'static str,
    /// Interpolated state at the event epoch.
    pub state: DynamicsState<S, C, F>,
}

/// Result returned by the propagation driver.
#[cfg(any(feature = "alloc", feature = "std"))]
#[derive(Debug, Clone)]
pub struct PropagationResult<S, C, F>
where
    S: ContinuousScale,
    C: ReferenceCenter,
    F: ReferenceFrame,
{
    /// Recorded samples, including the initial state and final state.
    pub samples: Vec<DynamicsState<S, C, F>>,
    /// Recorded event occurrences.
    pub events: Vec<EventOccurrence<S, C, F>>,
    /// Number of accepted steps.
    pub steps_taken: usize,
    /// Number of internally rejected trial steps reported by the adaptive integrator.
    pub steps_rejected: u32,
}

#[cfg(any(feature = "alloc", feature = "std"))]
fn lerp_state<S, C, F>(
    a: &DynamicsState<S, C, F>,
    b: &DynamicsState<S, C, F>,
    theta: f64,
) -> DynamicsState<S, C, F>
where
    S: ContinuousScale,
    C: ReferenceCenter,
    F: ReferenceFrame,
    C::Params: Clone,
{
    let px = a.position.x().value() + theta * (b.position.x().value() - a.position.x().value());
    let py = a.position.y().value() + theta * (b.position.y().value() - a.position.y().value());
    let pz = a.position.z().value() + theta * (b.position.z().value() - a.position.z().value());
    let vx = a.velocity.x().value() + theta * (b.velocity.x().value() - a.velocity.x().value());
    let vy = a.velocity.y().value() + theta * (b.velocity.y().value() - a.velocity.y().value());
    let vz = a.velocity.z().value() + theta * (b.velocity.z().value() - a.velocity.z().value());
    DynamicsState::new(
        a.epoch + Second::new(theta * (b.epoch - a.epoch).value()),
        affn::cartesian::Position::<C, F, qtty::unit::Kilometer>::new_with_params(
            a.position.center_params().clone(),
            px,
            py,
            pz,
        ),
        affn::cartesian::Velocity::<F, qtty::KmPerSecond>::new(vx, vy, vz),
    )
}

/// Propagate a state with an adaptive stepper over the interval configured in `cfg`.
#[cfg(any(feature = "alloc", feature = "std"))]
#[allow(clippy::too_many_lines)]
pub fn propagate<I, M, Ctx, S, C, F>(
    integrator: &I,
    model: &M,
    initial: DynamicsState<S, C, F>,
    cfg: &PropagationConfig<Ctx, S, C, F>,
    ctx: &Ctx,
) -> Result<PropagationResult<S, C, F>, PropagationError>
where
    I: AdaptiveStepper<Ctx, S, C, F>,
    M: AccelerationModel<Ctx, S, C, F>,
    S: ContinuousScale,
    C: ReferenceCenter,
    F: ReferenceFrame,
    C::Params: Clone,
{
    cfg.validate()
        .map_err(PropagationError::InvalidConfiguration)?;
    let total = cfg.total_duration_s();
    let direction = if total >= 0.0 { 1.0 } else { -1.0 };
    let mut samples = Vec::new();
    let mut events = Vec::new();
    let mut steps_taken = 0usize;
    let mut steps_rejected = 0u32;
    samples.push(initial.clone());
    if total == 0.0 {
        return Ok(PropagationResult {
            samples,
            events,
            steps_taken,
            steps_rejected,
        });
    }

    let mut g_prev = Vec::with_capacity(cfg.events.len());
    for ev in &cfg.events {
        g_prev.push(ev.evaluate(&initial, ctx).map_err(|source| {
            PropagationError::EventEvaluation {
                name: ev.name(),
                source,
            }
        })?);
    }

    let mut state = initial;
    let mut h = direction * cfg.h0.value().abs().min(cfg.h_max.value().abs());
    let mut next_output_t = cfg
        .output_every
        .map(|dt| cfg.t_start + Second::new(direction * dt.value().abs()));
    let mut output_at_iter = cfg.output_at.iter().peekable();
    let mut terminated = false;

    while !terminated {
        let remaining = (cfg.t_end - state.epoch).value();
        if remaining * direction <= 1e-12 {
            break;
        }
        if h.abs() > remaining.abs() {
            h = remaining;
        }
        if h.abs() > cfg.h_max.value().abs() {
            h = direction * cfg.h_max.value().abs();
        }
        if h.abs() < cfg.h_min.value().abs() && h.abs() < remaining.abs() {
            return Err(PropagationError::StepBelowMinimum {
                h_requested: h.abs(),
                h_min: cfg.h_min.value().abs(),
            });
        }

        if let Some(t_target) = next_output_t {
            let to_target = (t_target - state.epoch).value();
            if to_target.abs() > 1e-12 && to_target * direction >= 0.0 && to_target.abs() < h.abs()
            {
                h = to_target;
            }
        }
        if let Some(&t_target) = output_at_iter.peek() {
            let to_target = (*t_target - state.epoch).value();
            if to_target.abs() > 1e-12 && to_target * direction >= 0.0 && to_target.abs() < h.abs()
            {
                h = to_target;
            }
        }

        let (new_state, _h_used, h_next, rejected) = integrator
            .step(model, &state, Second::new(h), ctx)
            .map_err(PropagationError::StepControl)?;
        steps_taken += 1;
        steps_rejected += rejected;
        if steps_taken > cfg.max_steps {
            return Err(PropagationError::MaxStepsExceeded {
                max_steps: cfg.max_steps,
            });
        }

        for (i, ev) in cfg.events.iter().enumerate() {
            let g_new = ev.evaluate(&new_state, ctx).map_err(|source| {
                PropagationError::EventEvaluation {
                    name: ev.name(),
                    source,
                }
            })?;
            let sign_changed = g_prev[i] == 0.0 || g_new == 0.0 || (g_prev[i] * g_new) < 0.0;
            if sign_changed && ev.accepts_crossing(g_prev[i], g_new) {
                let theta = if (g_new - g_prev[i]).abs() > 1e-300 {
                    g_prev[i].abs() / (g_prev[i].abs() + g_new.abs())
                } else {
                    0.5
                };
                events.push(EventOccurrence {
                    event_name: ev.name(),
                    state: lerp_state(&state, &new_state, theta),
                });
                if ev.terminal() {
                    terminated = true;
                }
            }
            g_prev[i] = g_new;
        }

        if let Some(t_target) = next_output_t {
            let crossed = (t_target - state.epoch).value() * direction >= 0.0
                && (new_state.epoch - t_target).value() * direction >= 0.0;
            if crossed {
                samples.push(new_state.clone());
                if let Some(dt) = cfg.output_every {
                    next_output_t = Some(t_target + Second::new(direction * dt.value().abs()));
                }
            }
        }
        while let Some(&t_target) = output_at_iter.peek() {
            let crossed = (*t_target - state.epoch).value() * direction >= 0.0
                && (new_state.epoch - *t_target).value() * direction >= 0.0;
            if crossed {
                samples.push(new_state.clone());
                output_at_iter.next();
            } else {
                break;
            }
        }

        state = new_state;
        h = h_next.value();
    }

    if samples
        .last()
        .map(|sample| (sample.epoch - state.epoch).value().abs() > 0.0)
        .unwrap_or(true)
    {
        samples.push(state);
    }

    Ok(PropagationResult {
        samples,
        events,
        steps_taken,
        steps_rejected,
    })
}

#[cfg(all(test, any(feature = "alloc", feature = "std")))]
mod tests {
    use super::*;
    use crate::integrators::Dopri5;
    use crate::models::TwoBody;
    use affn::centers::ReferenceCenter;
    use affn::frames::ReferenceFrame;
    use qtty::{GravitationalParameter, Second};
    use tempoch::{Time, TT};

    #[derive(Debug, Clone, Copy)]
    struct Inertial;
    impl ReferenceFrame for Inertial {
        fn frame_name() -> &'static str {
            "Inertial"
        }
    }

    #[derive(Debug, Clone, Copy)]
    struct Center;
    impl ReferenceCenter for Center {
        type Params = ();
        fn center_name() -> &'static str {
            "Center"
        }
    }

    fn circular_state() -> DynamicsState<TT, Center, Inertial> {
        let mu = 398_600.441_8_f64;
        let r = 7000.0_f64;
        let v = (mu / r).sqrt();
        DynamicsState::new(
            Time::<TT>::from_raw_j2000_seconds(Second::new(0.0)).unwrap(),
            affn::cartesian::Position::<Center, Inertial, qtty::unit::Kilometer>::new(r, 0.0, 0.0),
            affn::cartesian::Velocity::<Inertial, qtty::KmPerSecond>::new(0.0, v, 0.0),
        )
    }

    #[test]
    fn validate_rejects_non_finite_span() {
        let epoch = Time::<TT>::from_raw_j2000_seconds(Second::new(0.0)).unwrap();
        let cfg = PropagationConfig::<(), TT, Center, Inertial>::new(
            epoch,
            epoch + Second::new(f64::NAN),
        );
        assert!(matches!(
            cfg.validate(),
            Err(PrincipiaError::InvalidPropagationConfig { .. })
        ));
    }

    #[test]
    fn zero_duration_returns_initial_state() {
        let s0 = circular_state();
        let cfg = PropagationConfig::<(), TT, Center, Inertial>::new(s0.epoch, s0.epoch);
        let result = propagate(
            &Dopri5::new(IntegratorTolerances::uniform(1e-9, 1e-6, 1e-9)),
            &TwoBody::new(GravitationalParameter::new(398_600.441_8)),
            s0,
            &cfg,
            &(),
        )
        .unwrap();
        assert_eq!(result.samples.len(), 1);
        assert_eq!(result.steps_taken, 0);
    }

    #[test]
    fn radial_threshold_event_fires() {
        let s0 = circular_state();
        let period = 2.0 * core::f64::consts::PI * (7000.0_f64.powi(3) / 398_600.441_8).sqrt();
        let ecc_state = DynamicsState::new(
            s0.epoch,
            s0.position,
            affn::cartesian::Velocity::<Inertial, qtty::KmPerSecond>::new(
                0.0,
                s0.velocity.y().value() * 1.01,
                0.0,
            ),
        );
        let cfg = PropagationConfig::<(), TT, Center, Inertial>::new(
            ecc_state.epoch,
            ecc_state.epoch + Second::new(period),
        )
        .with_event(Box::new(RadialThresholdEvent::new(Kilometers::new(7121.0))));
        let result = propagate(
            &Dopri5::new(IntegratorTolerances::uniform(1e-9, 1e-6, 1e-9)),
            &TwoBody::new(GravitationalParameter::new(398_600.441_8)),
            ecc_state,
            &cfg,
            &(),
        )
        .unwrap();
        assert!(!result.events.is_empty());
    }

    #[test]
    fn propagation_error_display() {
        let e = PropagationError::InvalidConfiguration(PrincipiaError::InvalidPropagationConfig {
            reason: "bad",
        });
        assert!(e.to_string().contains("bad"), "got: {e}");

        let e =
            PropagationError::StepControl(PrincipiaError::StepControlFailed { reason: "diverged" });
        assert!(e.to_string().contains("diverged"), "got: {e}");

        let e = PropagationError::StepBelowMinimum {
            h_requested: 1e-9,
            h_min: 1e-6,
        };
        assert!(
            e.to_string().contains("1e-9") || e.to_string().contains("1e"),
            "got: {e}"
        );

        let e = PropagationError::MaxStepsExceeded { max_steps: 42 };
        assert!(e.to_string().contains("42"), "got: {e}");

        let e = PropagationError::EventEvaluation {
            name: "my_event",
            source: PrincipiaError::DegenerateGeometry { reason: "r=0" },
        };
        assert!(e.to_string().contains("my_event"), "got: {e}");
    }

    #[test]
    fn propagation_error_source() {
        use core::error::Error;
        let e = PropagationError::InvalidConfiguration(PrincipiaError::InvalidPropagationConfig {
            reason: "bad",
        });
        assert!(e.source().is_some());

        let e = PropagationError::MaxStepsExceeded { max_steps: 1 };
        assert!(e.source().is_none());

        let e = PropagationError::StepBelowMinimum {
            h_requested: 1e-9,
            h_min: 1e-6,
        };
        assert!(e.source().is_none());
    }

    #[test]
    fn validate_rejects_invalid_h0() {
        let epoch = Time::<TT>::from_raw_j2000_seconds(Second::new(0.0)).unwrap();
        let mut cfg =
            PropagationConfig::<(), TT, Center, Inertial>::new(epoch, epoch + Second::new(100.0));
        cfg.h0 = Second::new(0.0);
        assert!(matches!(
            cfg.validate(),
            Err(PrincipiaError::InvalidPropagationConfig { .. })
        ));
    }

    #[test]
    fn validate_rejects_h_min_gt_h_max() {
        let epoch = Time::<TT>::from_raw_j2000_seconds(Second::new(0.0)).unwrap();
        let cfg =
            PropagationConfig::<(), TT, Center, Inertial>::new(epoch, epoch + Second::new(100.0))
                .with_min_step(Second::new(1000.0))
                .with_max_step(Second::new(10.0));
        assert!(matches!(
            cfg.validate(),
            Err(PrincipiaError::InvalidPropagationConfig { .. })
        ));
    }

    #[test]
    fn validate_rejects_max_steps_zero() {
        let epoch = Time::<TT>::from_raw_j2000_seconds(Second::new(0.0)).unwrap();
        let cfg =
            PropagationConfig::<(), TT, Center, Inertial>::new(epoch, epoch + Second::new(100.0))
                .with_max_steps(0);
        assert!(matches!(
            cfg.validate(),
            Err(PrincipiaError::InvalidPropagationConfig { .. })
        ));
    }

    #[test]
    fn validate_accepts_valid_config() {
        let epoch = Time::<TT>::from_raw_j2000_seconds(Second::new(0.0)).unwrap();
        let cfg =
            PropagationConfig::<(), TT, Center, Inertial>::new(epoch, epoch + Second::new(3600.0));
        assert!(cfg.validate().is_ok());
    }

    #[test]
    fn propagate_with_output_every_produces_samples() {
        let s0 = circular_state();
        let mu = 398_600.441_8;
        let r = 7000.0_f64;
        let period = 2.0 * core::f64::consts::PI * (r * r * r / mu).sqrt();
        let half = period / 2.0;
        let cfg = PropagationConfig::<(), TT, Center, Inertial>::new(
            s0.epoch,
            s0.epoch + Second::new(half),
        )
        .with_output_every(Second::new(half / 4.0))
        .with_initial_step(Second::new(60.0));
        let result = propagate(
            &Dopri5::new(IntegratorTolerances::uniform(1e-9, 1e-6, 1e-9)),
            &TwoBody::new(GravitationalParameter::new(mu)),
            s0,
            &cfg,
            &(),
        )
        .unwrap();
        // Should have initial + ~4 cadence samples + final
        assert!(
            result.samples.len() >= 3,
            "got {} samples",
            result.samples.len()
        );
    }

    #[test]
    fn propagate_with_output_at_explicit_epochs() {
        let s0 = circular_state();
        let mu = 398_600.441_8;
        let t1 = s0.epoch + Second::new(120.0);
        let t2 = s0.epoch + Second::new(240.0);
        let cfg = PropagationConfig::<(), TT, Center, Inertial>::new(
            s0.epoch,
            s0.epoch + Second::new(360.0),
        )
        .with_output_at(alloc::vec![t1, t2])
        .with_initial_step(Second::new(60.0));
        let result = propagate(
            &Dopri5::new(IntegratorTolerances::uniform(1e-9, 1e-6, 1e-9)),
            &TwoBody::new(GravitationalParameter::new(mu)),
            s0,
            &cfg,
            &(),
        )
        .unwrap();
        assert!(result.samples.len() >= 3);
    }

    #[test]
    fn propagate_terminal_event_stops_early() {
        let s0 = circular_state();
        let mu = 398_600.441_8;
        let r = 7000.0_f64;
        let period = 2.0 * core::f64::consts::PI * (r * r * r / mu).sqrt();
        // Slightly eccentric orbit crosses 7050 km on the way up
        let ecc_state = DynamicsState::new(
            s0.epoch,
            s0.position,
            affn::cartesian::Velocity::<Inertial, qtty::KmPerSecond>::new(
                0.0,
                s0.velocity.y().value() * 1.005,
                0.0,
            ),
        );
        let event = RadialThresholdEvent::new(Kilometers::new(7050.0))
            .terminal(true)
            .direction(EventDirection::Rising);
        let cfg = PropagationConfig::<(), TT, Center, Inertial>::new(
            ecc_state.epoch,
            ecc_state.epoch + Second::new(period),
        )
        .with_event(Box::new(event));
        let result = propagate(
            &Dopri5::new(IntegratorTolerances::uniform(1e-9, 1e-6, 1e-9)),
            &TwoBody::new(GravitationalParameter::new(mu)),
            ecc_state,
            &cfg,
            &(),
        )
        .unwrap();
        // Terminal event should have stopped propagation before one full period
        assert!(
            result.steps_taken < 10_000,
            "steps={} — expected early termination",
            result.steps_taken
        );
        assert!(!result.events.is_empty());
    }

    #[test]
    fn max_steps_exceeded_returns_error() {
        let s0 = circular_state();
        let mu = 398_600.441_8;
        let r = 7000.0_f64;
        let period = 2.0 * core::f64::consts::PI * (r * r * r / mu).sqrt();
        let cfg = PropagationConfig::<(), TT, Center, Inertial>::new(
            s0.epoch,
            s0.epoch + Second::new(period),
        )
        .with_max_steps(1); // allow only 1 step
        let result = propagate(
            &Dopri5::new(IntegratorTolerances::uniform(1e-9, 1e-6, 1e-9)),
            &TwoBody::new(GravitationalParameter::new(mu)),
            s0,
            &cfg,
            &(),
        );
        assert!(
            matches!(result, Err(PropagationError::MaxStepsExceeded { .. })),
            "expected MaxStepsExceeded, got: {result:?}"
        );
    }

    #[test]
    fn event_direction_rising_ignores_falling() {
        // State starts above threshold 7000 and orbital dynamics keep it there in one step.
        // A falling event (from above) should not fire.
        let s0 = circular_state();
        let mu = 398_600.441_8;
        // Use a threshold below the orbit (6800 km) with Rising direction → shouldn't fire
        let event =
            RadialThresholdEvent::new(Kilometers::new(6800.0)).direction(EventDirection::Rising);
        let cfg = PropagationConfig::<(), TT, Center, Inertial>::new(
            s0.epoch,
            s0.epoch + Second::new(600.0),
        )
        .with_event(Box::new(event));
        let result = propagate(
            &Dopri5::new(IntegratorTolerances::uniform(1e-9, 1e-6, 1e-9)),
            &TwoBody::new(GravitationalParameter::new(mu)),
            s0,
            &cfg,
            &(),
        )
        .unwrap();
        // No rising crossing of 6800 km from a 7000 km circular orbit
        assert!(
            result.events.is_empty(),
            "no rising crossing expected, got {}",
            result.events.len()
        );
    }

    #[test]
    fn event_direction_either_fires_on_any_crossing() {
        let s0 = circular_state();
        let mu = 398_600.441_8;
        let r = 7000.0_f64;
        let period = 2.0 * core::f64::consts::PI * (r * r * r / mu).sqrt();
        // Slightly eccentric orbit — will cross a threshold both ways
        let ecc_state = DynamicsState::new(
            s0.epoch,
            s0.position,
            affn::cartesian::Velocity::<Inertial, qtty::KmPerSecond>::new(
                0.0,
                s0.velocity.y().value() * 1.01,
                0.0,
            ),
        );
        let event =
            RadialThresholdEvent::new(Kilometers::new(7120.0)).direction(EventDirection::Either);
        let cfg = PropagationConfig::<(), TT, Center, Inertial>::new(
            ecc_state.epoch,
            ecc_state.epoch + Second::new(period),
        )
        .with_event(Box::new(event));
        let result = propagate(
            &Dopri5::new(IntegratorTolerances::uniform(1e-9, 1e-6, 1e-9)),
            &TwoBody::new(GravitationalParameter::new(mu)),
            ecc_state,
            &cfg,
            &(),
        )
        .unwrap();
        assert!(!result.events.is_empty());
    }

    #[test]
    fn propagate_builder_methods_work() {
        let epoch = Time::<TT>::from_raw_j2000_seconds(Second::new(0.0)).unwrap();
        let cfg =
            PropagationConfig::<(), TT, Center, Inertial>::new(epoch, epoch + Second::new(100.0))
                .with_initial_step(Second::new(5.0))
                .with_max_step(Second::new(50.0))
                .with_min_step(Second::new(0.5))
                .with_tolerances(IntegratorTolerances::uniform(1e-8, 1e-5, 1e-8))
                .with_max_steps(500);
        assert_eq!(cfg.h0.value(), 5.0);
        assert_eq!(cfg.h_max.value(), 50.0);
        assert_eq!(cfg.h_min.value(), 0.5);
        assert_eq!(cfg.max_steps, 500);
        assert!((cfg.total_duration_s() - 100.0).abs() < 1e-12);
    }

    #[test]
    fn radial_threshold_event_name() {
        let ev = RadialThresholdEvent::new(Kilometers::new(7000.0));
        let state = circular_state();
        assert_eq!(
            <RadialThresholdEvent as EventDetector<(), TT, Center, Inertial>>::name(&ev),
            "radial_threshold"
        );
        let g = <RadialThresholdEvent as EventDetector<(), TT, Center, Inertial>>::evaluate(
            &ev,
            &state,
            &(),
        )
        .unwrap();
        // Orbit is at exactly 7000 km, threshold is 7000 km → g ≈ 0
        assert!(g.abs() < 1e-6, "g={g}");
    }

    #[test]
    fn event_direction_falling_fires_correctly() {
        let s0 = circular_state();
        let mu = 398_600.441_8;
        let r = 7000.0_f64;
        let period = 2.0 * core::f64::consts::PI * (r * r * r / mu).sqrt();
        let ecc_state = DynamicsState::new(
            s0.epoch,
            s0.position,
            affn::cartesian::Velocity::<Inertial, qtty::KmPerSecond>::new(
                0.0,
                s0.velocity.y().value() * 1.01,
                0.0,
            ),
        );
        // Falling: orbit approaches from above apogee
        let event =
            RadialThresholdEvent::new(Kilometers::new(7121.0)).direction(EventDirection::Falling);
        let cfg = PropagationConfig::<(), TT, Center, Inertial>::new(
            ecc_state.epoch,
            ecc_state.epoch + Second::new(period),
        )
        .with_event(Box::new(event));
        let result = propagate(
            &Dopri5::new(IntegratorTolerances::uniform(1e-9, 1e-6, 1e-9)),
            &TwoBody::new(GravitationalParameter::new(mu)),
            ecc_state,
            &cfg,
            &(),
        )
        .unwrap();
        // Should see a falling crossing on the way back down
        assert!(!result.events.is_empty());
    }
}