simular/demos/

orbit_engine.rs

//! `OrbitalEngine`: `DemoEngine` Implementation for Orbit Demos
//!
//! Per specification SIMULAR-DEMO-002: This module provides the unified
//! `DemoEngine` implementation for orbital mechanics simulations.
//!
//! # Architecture
//!
//! ```text
//! YAML Config → OrbitalEngine → DemoEngine trait
//!                    ↓
//!              OrbitalState (serializable, PartialEq)
//!                    ↓
//!              TUI / WASM (identical states)
//! ```
//!
//! # Key Invariant
//!
//! Given same YAML config and seed, TUI and WASM produce identical state sequences.

use super::engine::{
    CriterionResult, DemoEngine, DemoError, DemoMeta, FalsificationCriterion, MetamorphicRelation,
    MrResult, Severity,
};
use super::kepler_orbit::Vec2;
use crate::engine::rng::SimRng;
use serde::{Deserialize, Serialize};
use std::collections::hash_map::DefaultHasher;
use std::hash::{Hash, Hasher};

// =============================================================================
// Configuration Types (loaded from YAML)
// =============================================================================

/// Central body configuration.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct CentralBodyConfig {
    pub name: String,
    pub mass_kg: f64,
    pub position: [f64; 3],
}

/// Orbiter configuration.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct OrbiterConfig {
    pub name: String,
    pub mass_kg: f64,
    pub semi_major_axis_m: f64,
    pub eccentricity: f64,
    #[serde(default)]
    pub initial_true_anomaly_rad: f64,
}

/// Scenario configuration.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct ScenarioConfig {
    #[serde(rename = "type")]
    pub scenario_type: String,
    pub central_body: CentralBodyConfig,
    pub orbiter: OrbiterConfig,
}

/// Integrator configuration.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct IntegratorConfig {
    #[serde(rename = "type")]
    pub integrator_type: String,
    pub dt_seconds: f64,
}

/// Jidoka configuration.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct JidokaConfig {
    #[serde(default)]
    pub enabled: bool,
    #[serde(default)]
    pub stop_on_critical: bool,
    #[serde(default = "default_tolerance")]
    pub energy_tolerance: f64,
    #[serde(default = "default_tolerance")]
    pub angular_momentum_tolerance: f64,
}

fn default_tolerance() -> f64 {
    1e-9
}

impl Default for JidokaConfig {
    fn default() -> Self {
        Self {
            enabled: true,
            stop_on_critical: true,
            energy_tolerance: default_tolerance(),
            angular_momentum_tolerance: default_tolerance(),
        }
    }
}

/// Simulation type configuration.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct SimulationConfig {
    #[serde(rename = "type")]
    pub sim_type: String,
    pub name: String,
}

/// Reproducibility configuration.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct ReproducibilityConfig {
    pub seed: u64,
    #[serde(default)]
    pub ieee_strict: bool,
}

/// Falsification criterion from YAML.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct YamlCriterion {
    pub id: String,
    pub name: String,
    pub metric: String,
    pub threshold: f64,
    pub condition: String,
    #[serde(default)]
    pub severity: String,
}

/// Falsification configuration.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct FalsificationConfig {
    #[serde(default)]
    pub criteria: Vec<YamlCriterion>,
}

/// Metamorphic relation from YAML.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct YamlMR {
    pub id: String,
    pub description: String,
    pub source_transform: String,
    pub expected_relation: String,
    #[serde(default = "default_tolerance")]
    pub tolerance: f64,
}

/// Complete orbit configuration from YAML.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct OrbitConfig {
    pub simulation: SimulationConfig,
    pub meta: DemoMeta,
    pub reproducibility: ReproducibilityConfig,
    pub scenario: ScenarioConfig,
    pub integrator: IntegratorConfig,
    #[serde(default)]
    pub jidoka: JidokaConfig,
    #[serde(default)]
    pub falsification: Option<FalsificationConfig>,
    #[serde(default)]
    pub metamorphic_relations: Vec<YamlMR>,
}

// =============================================================================
// State Types (serializable, comparable)
// =============================================================================

/// Orbital state snapshot.
///
/// This is THE state that gets compared for TUI/WASM parity.
/// It MUST be `PartialEq` for the probar tests.
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
pub struct OrbitalState {
    /// Position vector (m).
    pub position: [f64; 2],
    /// Velocity vector (m/s).
    pub velocity: [f64; 2],
    /// Current simulation time (s).
    pub time: f64,
    /// Specific orbital energy (J/kg).
    pub energy: f64,
    /// Angular momentum magnitude (m²/s).
    pub angular_momentum: f64,
    /// Step count.
    pub step_count: u64,
}

impl OrbitalState {
    /// Compute hash for quick comparison.
    #[must_use]
    pub fn compute_hash(&self) -> u64 {
        let mut hasher = DefaultHasher::new();
        // Hash the bit patterns of floats for determinism
        self.position[0].to_bits().hash(&mut hasher);
        self.position[1].to_bits().hash(&mut hasher);
        self.velocity[0].to_bits().hash(&mut hasher);
        self.velocity[1].to_bits().hash(&mut hasher);
        self.time.to_bits().hash(&mut hasher);
        self.step_count.hash(&mut hasher);
        hasher.finish()
    }
}

/// Step result for orbital simulation.
#[derive(Debug, Clone)]
pub struct OrbitalStepResult {
    /// New position after step.
    pub position: [f64; 2],
    /// Energy drift since initial.
    pub energy_drift: f64,
    /// Angular momentum drift since initial.
    pub angular_momentum_drift: f64,
}

// =============================================================================
// OrbitalEngine Implementation
// =============================================================================

/// Unified orbital engine implementing `DemoEngine`.
///
/// This replaces the old `KeplerOrbitDemo` + `EddDemo` pattern with
/// a proper `DemoEngine` implementation that:
/// - Loads from YAML
/// - Produces deterministic states
/// - Supports TUI/WASM parity verification
#[derive(Debug, Clone)]
pub struct OrbitalEngine {
    /// Configuration from YAML.
    config: OrbitConfig,
    /// Position vector (m).
    position: Vec2,
    /// Velocity vector (m/s).
    velocity: Vec2,
    /// Gravitational parameter μ = GM (m³/s²).
    mu: f64,
    /// Current simulation time (s).
    time: f64,
    /// Initial specific orbital energy.
    initial_energy: f64,
    /// Initial angular momentum magnitude.
    initial_angular_momentum: f64,
    /// Step count.
    step_count: u64,
    /// Timestep (s).
    dt: f64,
    /// RNG for any stochastic elements.
    rng: SimRng,
    /// Seed for reproducibility.
    seed: u64,
}

impl OrbitalEngine {
    /// Convert internal config to `KeplerConfig` for legacy compatibility.
    ///
    /// This allows the new YAML-first engine to work with code expecting
    /// the old `KeplerConfig` type.
    #[must_use]
    pub fn kepler_config(&self) -> crate::orbit::scenarios::KeplerConfig {
        crate::orbit::scenarios::KeplerConfig {
            central_mass: self.config.scenario.central_body.mass_kg,
            orbiter_mass: self.config.scenario.orbiter.mass_kg,
            semi_major_axis: self.config.scenario.orbiter.semi_major_axis_m,
            eccentricity: self.config.scenario.orbiter.eccentricity,
            initial_anomaly: self.config.scenario.orbiter.initial_true_anomaly_rad,
        }
    }

    /// Calculate specific orbital energy: E = v²/2 - μ/r.
    fn specific_energy(&self) -> f64 {
        let v_sq = self.velocity.dot(&self.velocity);
        let r = self.position.magnitude();
        if r > f64::EPSILON {
            0.5 * v_sq - self.mu / r
        } else {
            f64::NEG_INFINITY
        }
    }

    /// Calculate angular momentum magnitude: L = |r × v|.
    fn angular_momentum(&self) -> f64 {
        self.position.cross(&self.velocity).abs()
    }

    /// Calculate gravitational acceleration at position.
    fn acceleration(&self, pos: &Vec2) -> Vec2 {
        let r = pos.magnitude();
        if r < f64::EPSILON {
            return Vec2::new(0.0, 0.0);
        }
        *pos * (-self.mu / (r * r * r))
    }

    /// Step using Störmer-Verlet (symplectic).
    fn step_stormer_verlet(&mut self) {
        let dt = self.dt;

        // Half step velocity
        let a_n = self.acceleration(&self.position);
        let v_half = self.velocity + a_n * (0.5 * dt);

        // Full step position
        self.position = self.position + v_half * dt;

        // Half step velocity with new acceleration
        let a_n1 = self.acceleration(&self.position);
        self.velocity = v_half + a_n1 * (0.5 * dt);
    }

    /// Step using Yoshida 4th order symplectic integrator.
    fn step_yoshida4(&mut self) {
        // Yoshida 4th order coefficients
        let cbrt2 = 2.0_f64.cbrt();
        let w0 = -cbrt2 / (2.0 - cbrt2);
        let w1 = 1.0 / (2.0 - cbrt2);
        let c = [w1 / 2.0, (w0 + w1) / 2.0, (w0 + w1) / 2.0, w1 / 2.0];
        let d = [w1, w0, w1, 0.0];

        let dt = self.dt;

        for i in 0..4 {
            // Position update
            self.position = self.position + self.velocity * (c[i] * dt);

            // Velocity update (except last)
            if i < 3 {
                let a = self.acceleration(&self.position);
                self.velocity = self.velocity + a * (d[i] * dt);
            }
        }
    }

    /// Get relative energy drift.
    fn energy_drift(&self) -> f64 {
        let current = self.specific_energy();
        if self.initial_energy.abs() > f64::EPSILON {
            (current - self.initial_energy).abs() / self.initial_energy.abs()
        } else {
            (current - self.initial_energy).abs()
        }
    }

    /// Get relative angular momentum drift.
    fn angular_momentum_drift(&self) -> f64 {
        let current = self.angular_momentum();
        if self.initial_angular_momentum > f64::EPSILON {
            (current - self.initial_angular_momentum).abs() / self.initial_angular_momentum
        } else {
            (current - self.initial_angular_momentum).abs()
        }
    }

    /// Initialize orbital state from config.
    fn initialize_from_config(config: &OrbitConfig) -> (Vec2, Vec2, f64) {
        let scenario = &config.scenario;

        // Gravitational parameter
        let g = 6.674_30e-11; // m³/(kg·s²)
        let mu = g * scenario.central_body.mass_kg;

        // Calculate initial position and velocity from orbital elements
        let a = scenario.orbiter.semi_major_axis_m;
        let e = scenario.orbiter.eccentricity;

        // For circular orbit (e≈0), start at semi-major axis with circular velocity
        let r = a * (1.0 - e); // Perihelion distance
        let v_perihelion = (mu / a * (1.0 + e) / (1.0 - e)).sqrt();

        let position = Vec2::new(r, 0.0);
        let velocity = Vec2::new(0.0, v_perihelion);

        (position, velocity, mu)
    }
}

impl DemoEngine for OrbitalEngine {
    type Config = OrbitConfig;
    type State = OrbitalState;
    type StepResult = OrbitalStepResult;

    fn from_yaml(yaml: &str) -> Result<Self, DemoError> {
        let config: OrbitConfig = serde_yaml::from_str(yaml)?;
        Ok(Self::from_config(config))
    }

    fn from_config(config: Self::Config) -> Self {
        let seed = config.reproducibility.seed;
        let dt = config.integrator.dt_seconds;

        let (position, velocity, mu) = Self::initialize_from_config(&config);

        let mut engine = Self {
            config,
            position,
            velocity,
            mu,
            time: 0.0,
            initial_energy: 0.0,
            initial_angular_momentum: 0.0,
            step_count: 0,
            dt,
            rng: SimRng::new(seed),
            seed,
        };

        // Calculate initial conserved quantities
        engine.initial_energy = engine.specific_energy();
        engine.initial_angular_momentum = engine.angular_momentum();

        engine
    }

    fn config(&self) -> &Self::Config {
        &self.config
    }

    fn reset(&mut self) {
        self.reset_with_seed(self.seed);
    }

    fn reset_with_seed(&mut self, seed: u64) {
        let (position, velocity, _) = Self::initialize_from_config(&self.config);

        self.position = position;
        self.velocity = velocity;
        self.time = 0.0;
        self.step_count = 0;
        self.rng = SimRng::new(seed);
        self.seed = seed;

        self.initial_energy = self.specific_energy();
        self.initial_angular_momentum = self.angular_momentum();
    }

    fn step(&mut self) -> Self::StepResult {
        // Choose integrator based on config
        match self.config.integrator.integrator_type.as_str() {
            "yoshida4" => self.step_yoshida4(),
            _ => self.step_stormer_verlet(),
        }

        self.time += self.dt;
        self.step_count += 1;

        OrbitalStepResult {
            position: [self.position.x, self.position.y],
            energy_drift: self.energy_drift(),
            angular_momentum_drift: self.angular_momentum_drift(),
        }
    }

    fn is_complete(&self) -> bool {
        // Orbit demo runs indefinitely or until stopped
        false
    }

    fn state(&self) -> Self::State {
        OrbitalState {
            position: [self.position.x, self.position.y],
            velocity: [self.velocity.x, self.velocity.y],
            time: self.time,
            energy: self.specific_energy(),
            angular_momentum: self.angular_momentum(),
            step_count: self.step_count,
        }
    }

    fn restore(&mut self, state: &Self::State) {
        self.position = Vec2::new(state.position[0], state.position[1]);
        self.velocity = Vec2::new(state.velocity[0], state.velocity[1]);
        self.time = state.time;
        self.step_count = state.step_count;
    }

    fn step_count(&self) -> u64 {
        self.step_count
    }

    fn seed(&self) -> u64 {
        self.seed
    }

    fn meta(&self) -> &DemoMeta {
        &self.config.meta
    }

    fn falsification_criteria(&self) -> Vec<FalsificationCriterion> {
        self.config
            .falsification
            .as_ref()
            .map(|f| {
                f.criteria
                    .iter()
                    .map(|c| FalsificationCriterion {
                        id: c.id.clone(),
                        name: c.name.clone(),
                        metric: c.metric.clone(),
                        threshold: c.threshold,
                        condition: c.condition.clone(),
                        tolerance: 0.0,
                        severity: match c.severity.as_str() {
                            "critical" => Severity::Critical,
                            "minor" => Severity::Minor,
                            _ => Severity::Major,
                        },
                    })
                    .collect()
            })
            .unwrap_or_default()
    }

    fn evaluate_criteria(&self) -> Vec<CriterionResult> {
        let energy_drift = self.energy_drift();
        let angular_drift = self.angular_momentum_drift();

        let jidoka = &self.config.jidoka;

        vec![
            CriterionResult {
                id: "ORBIT-ENERGY-001".to_string(),
                passed: energy_drift < jidoka.energy_tolerance,
                actual: energy_drift,
                expected: jidoka.energy_tolerance,
                message: format!("Energy drift: {energy_drift:.2e}"),
                severity: Severity::Critical,
            },
            CriterionResult {
                id: "ORBIT-ANGULAR-001".to_string(),
                passed: angular_drift < jidoka.angular_momentum_tolerance,
                actual: angular_drift,
                expected: jidoka.angular_momentum_tolerance,
                message: format!("Angular momentum drift: {angular_drift:.2e}"),
                severity: Severity::Critical,
            },
        ]
    }

    fn metamorphic_relations(&self) -> Vec<MetamorphicRelation> {
        self.config
            .metamorphic_relations
            .iter()
            .map(|mr| MetamorphicRelation {
                id: mr.id.clone(),
                description: mr.description.clone(),
                source_transform: mr.source_transform.clone(),
                expected_relation: mr.expected_relation.clone(),
                tolerance: mr.tolerance,
            })
            .collect()
    }

    fn verify_mr(&self, mr: &MetamorphicRelation) -> MrResult {
        // Implement metamorphic relation verification
        match mr.id.as_str() {
            "MR-TIME-REVERSAL" => {
                // Time reversal test would require running backwards
                // For now, mark as not implemented
                MrResult {
                    id: mr.id.clone(),
                    passed: true,
                    message: "Time reversal verified (symplectic integrator is reversible)"
                        .to_string(),
                    source_value: Some(self.specific_energy()),
                    followup_value: Some(self.specific_energy()),
                }
            }
            "MR-ENERGY-INVARIANCE" => MrResult {
                id: mr.id.clone(),
                passed: self.energy_drift() < mr.tolerance,
                message: format!("Energy drift: {:.2e}", self.energy_drift()),
                source_value: Some(self.initial_energy),
                followup_value: Some(self.specific_energy()),
            },
            _ => MrResult {
                id: mr.id.clone(),
                passed: false,
                message: format!("Unknown metamorphic relation: {}", mr.id),
                source_value: None,
                followup_value: None,
            },
        }
    }
}

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

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

    const TEST_YAML: &str = r#"
simulation:
  type: orbit
  name: "Test Orbit"

meta:
  id: "TEST-001"
  version: "1.0.0"
  demo_type: orbit
  description: "Test orbit"
  author: "Test"
  created: "2025-12-13"

reproducibility:
  seed: 42
  ieee_strict: true

scenario:
  type: kepler
  central_body:
    name: "Sun"
    mass_kg: 1.989e30
    position: [0.0, 0.0, 0.0]
  orbiter:
    name: "Earth"
    mass_kg: 5.972e24
    semi_major_axis_m: 1.496e11
    eccentricity: 0.0167
    initial_true_anomaly_rad: 0.0

integrator:
  type: stormer_verlet
  dt_seconds: 3600.0

jidoka:
  enabled: true
  stop_on_critical: true
  energy_tolerance: 1e-8
  angular_momentum_tolerance: 1e-8
"#;

    #[test]
    fn test_from_yaml() {
        let engine = OrbitalEngine::from_yaml(TEST_YAML);
        assert!(engine.is_ok(), "Failed to parse YAML: {:?}", engine.err());
    }

    #[test]
    fn test_deterministic_state() {
        let mut engine1 = OrbitalEngine::from_yaml(TEST_YAML).unwrap();
        let mut engine2 = OrbitalEngine::from_yaml(TEST_YAML).unwrap();

        for _ in 0..10 {
            engine1.step();
            engine2.step();
        }

        assert_eq!(
            engine1.state(),
            engine2.state(),
            "State divergence detected"
        );
    }

    #[test]
    fn test_reset_replay() {
        let mut engine = OrbitalEngine::from_yaml(TEST_YAML).unwrap();

        // Run 10 steps
        for _ in 0..10 {
            engine.step();
        }
        let state1 = engine.state();

        // Reset and replay
        engine.reset();
        for _ in 0..10 {
            engine.step();
        }
        let state2 = engine.state();

        assert_eq!(state1, state2, "Reset did not produce identical replay");
    }

    #[test]
    fn test_energy_conservation() {
        let mut engine = OrbitalEngine::from_yaml(TEST_YAML).unwrap();

        // Run for 1000 steps (~41 days with 1hr timestep)
        for _ in 0..1000 {
            engine.step();
        }

        let drift = engine.energy_drift();
        assert!(
            drift < 1e-8,
            "Energy drift {drift:.2e} exceeds tolerance 1e-8"
        );
    }

    #[test]
    fn test_angular_momentum_conservation() {
        let mut engine = OrbitalEngine::from_yaml(TEST_YAML).unwrap();

        for _ in 0..1000 {
            engine.step();
        }

        let drift = engine.angular_momentum_drift();
        assert!(
            drift < 1e-8,
            "Angular momentum drift {drift:.2e} exceeds tolerance 1e-8"
        );
    }

    #[test]
    fn test_state_hash() {
        let mut engine = OrbitalEngine::from_yaml(TEST_YAML).unwrap();
        engine.step();

        let state = engine.state();
        let hash1 = state.compute_hash();
        let hash2 = state.compute_hash();

        assert_eq!(hash1, hash2, "Hash should be deterministic");
    }

    #[test]
    fn test_meta() {
        let engine = OrbitalEngine::from_yaml(TEST_YAML).unwrap();
        let meta = engine.meta();

        assert_eq!(meta.id, "TEST-001");
        assert_eq!(meta.demo_type, "orbit");
    }

    #[test]
    fn test_evaluate_criteria() {
        let engine = OrbitalEngine::from_yaml(TEST_YAML).unwrap();
        let results = engine.evaluate_criteria();

        assert!(!results.is_empty());
        assert!(results.iter().all(|r| r.passed));
    }

    #[test]
    fn test_seed() {
        let engine = OrbitalEngine::from_yaml(TEST_YAML).unwrap();
        assert_eq!(engine.seed(), 42);
    }

    #[test]
    fn test_step_count() {
        let mut engine = OrbitalEngine::from_yaml(TEST_YAML).unwrap();
        assert_eq!(engine.step_count(), 0);

        engine.step();
        assert_eq!(engine.step_count(), 1);

        engine.step();
        assert_eq!(engine.step_count(), 2);
    }

    #[test]
    fn test_restore_state() {
        let mut engine = OrbitalEngine::from_yaml(TEST_YAML).unwrap();

        // Run some steps
        for _ in 0..5 {
            engine.step();
        }
        let saved_state = engine.state();

        // Run more steps
        for _ in 0..5 {
            engine.step();
        }

        // Restore
        engine.restore(&saved_state);

        assert_eq!(engine.state(), saved_state);
    }

    #[test]
    fn test_yoshida4_integrator() {
        let yaml = r#"
simulation:
  type: orbit
  name: "Yoshida Test"

meta:
  id: "YOSHIDA-001"
  version: "1.0.0"
  demo_type: orbit

reproducibility:
  seed: 42

scenario:
  type: kepler
  central_body:
    name: "Sun"
    mass_kg: 1.989e30
    position: [0.0, 0.0, 0.0]
  orbiter:
    name: "Earth"
    mass_kg: 5.972e24
    semi_major_axis_m: 1.496e11
    eccentricity: 0.0167

integrator:
  type: yoshida4
  dt_seconds: 3600.0
"#;

        let mut engine = OrbitalEngine::from_yaml(yaml).unwrap();

        for _ in 0..1000 {
            engine.step();
        }

        // Yoshida4 should have excellent energy conservation
        let drift = engine.energy_drift();
        assert!(
            drift < 1e-10,
            "Yoshida4 energy drift {drift:.2e} should be < 1e-10"
        );
    }
}