kshana 0.37.0

// SPDX-License-Identifier: AGPL-3.0-only
//! Impulsive / finite-burn maneuvers, a Lambert two-body transfer solver, and a
//! porkchop (departure × arrival C3 / arrival-V∞) sweep — the maneuver-modeling and
//! trajectory-design beachhead.
//!
//! Scope (honest). This is the *performance-simulation* layer above a full mission-design
//! tool (GMAT / Orekit), not a replacement for one. It ships:
//!   * impulsive ΔV nodes that apply a velocity discontinuity and carry a 6×6 covariance
//!     forward (deterministic burn ⇒ identity state-transition across the instant; the
//!     stochastic execution-error covariance adds into the velocity block, rotated from the
//!     burn frame),
//!   * a finite-burn integration (constant thrust over a burn arc, mass as a state) whose
//!     achieved ΔV is checked against the closed-form **Tsiolkovsky** rocket equation,
//!   * an **Izzo-2015** single-revolution Lambert solver (`r1`, `r2`, time-of-flight ⇒
//!     `v1`, `v2` for a two-body transfer), and
//!   * a **porkchop** sweep that maps a launch-epoch × arrival-epoch grid to departure C3
//!     and arrival V∞, emitted as a 2-D JSON array for browser contour rendering.
//!
//! Validation is fully self-contained and closed-form, **stronger** than reading a value off
//! a GMAT tutorial: every Lambert output is round-tripped through an exact universal-variable
//! Kepler propagator (two-body truth to machine precision — it must land back on `r2`), the
//! finite burn is checked against Tsiolkovsky to < 0.01 %, and the porkchop minimum is checked
//! against the analytic Hohmann-transfer C3 floor for two coplanar circular orbits.
//!
//! Honest residuals (not claimed): multi-revolution (M ≥ 1) Lambert branches, a planetary
//! ephemeris (the porkchop uses a synthetic coplanar-circular heliocentric model so the
//! optimum is analytic — a real DE-ephemeris Earth–Mars C3 cross-check against GMAT needs the
//! same ephemeris and has not been run), and live wiring of the porkchop JSON into the web
//! playground contour widget. Kshana points users to GMAT/Orekit for full multi-burn /
//! low-thrust optimization.

use crate::portable_math::PortableFloat;
use serde::Serialize;
use std::f64::consts::PI;

type Vec3 = [f64; 3];

/// Standard gravity used to convert specific impulse (s) to effective exhaust velocity (m/s).
pub const G0: f64 = 9.806_65;
/// Heliocentric gravitational parameter (m³/s²), IAU 2015 nominal `GM_Sun`.
pub const MU_SUN: f64 = 1.327_124_400_18e20;
/// Astronomical unit (m), IAU 2012 definition.
pub const AU_M: f64 = 1.495_978_707e11;

// ---------------------------------------------------------------------------
// Small vector helpers (kept local, mirroring the other modules' style).
// ---------------------------------------------------------------------------
fn add(a: Vec3, b: Vec3) -> Vec3 {
    [a[0] + b[0], a[1] + b[1], a[2] + b[2]]
}
fn sub(a: Vec3, b: Vec3) -> Vec3 {
    [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
}
fn scale(a: Vec3, s: f64) -> Vec3 {
    [a[0] * s, a[1] * s, a[2] * s]
}
fn dot(a: Vec3, b: Vec3) -> f64 {
    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
}
fn cross(a: Vec3, b: Vec3) -> Vec3 {
    [
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0],
    ]
}
fn norm(a: Vec3) -> f64 {
    dot(a, a).sqrt()
}
fn unit(a: Vec3) -> Vec3 {
    let n = norm(a);
    [a[0] / n, a[1] / n, a[2] / n]
}

// ===========================================================================
// 1. Impulsive maneuver with covariance propagation.
// ===========================================================================

/// The frame an impulsive ΔV (and its execution-error covariance) is expressed in.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum ManeuverFrame {
    /// Inertial: ΔV components are already in the propagation (ECI) frame.
    Eci,
    /// Local-vertical/local-horizontal (radial, along-track, cross-track) at the burn state.
    Lvlh,
}

/// An instantaneous velocity change with an optional execution-error covariance.
#[derive(Clone, Copy, Debug)]
pub struct ImpulsiveManeuver {
    /// ΔV components (m/s) expressed in [`ImpulsiveManeuver::frame`].
    pub dv: Vec3,
    /// The frame `dv` and `exec_cov` are given in.
    pub frame: ManeuverFrame,
    /// 3×3 ΔV execution-error covariance (m²/s²) in the same frame (zeros ⇒ perfect burn).
    pub exec_cov: [[f64; 3]; 3],
}

/// Columns of the LVLH→ECI rotation `[r̂ | t̂ | ĥ]` at state `(r, v)`:
/// radial `r̂ = r/|r|`, cross-track `ĥ = (r×v)/|r×v|`, along-track `t̂ = ĥ×r̂` (toward `+v`).
pub fn lvlh_to_eci(r: Vec3, v: Vec3) -> [[f64; 3]; 3] {
    let ur = unit(r);
    let uh = unit(cross(r, v));
    let ut = cross(uh, ur);
    // Column-major: m[i] is the i-th basis vector expressed in ECI.
    [ur, ut, uh]
}

/// Rotate a body-frame vector `x` into ECI given the LVLH→ECI columns `m`.
fn rotate(m: &[[f64; 3]; 3], x: Vec3) -> Vec3 {
    [
        m[0][0] * x[0] + m[1][0] * x[1] + m[2][0] * x[2],
        m[0][1] * x[0] + m[1][1] * x[1] + m[2][1] * x[2],
        m[0][2] * x[0] + m[1][2] * x[1] + m[2][2] * x[2],
    ]
}

/// Apply an impulsive maneuver to a 6-state `[r(3), v(3)]` and its 6×6 covariance.
///
/// Position is continuous; velocity jumps by the (frame-resolved) ΔV. Because the burn is a
/// deterministic, state-independent shift, the state-transition across the instant is the
/// identity, so the *a priori* covariance is unchanged — only the execution-error covariance
/// (rotated into ECI) adds into the velocity–velocity block. Returns `(state_after, cov_after)`.
pub fn apply_impulse(
    state: [f64; 6],
    cov: [[f64; 6]; 6],
    man: &ImpulsiveManeuver,
) -> ([f64; 6], [[f64; 6]; 6]) {
    let r = [state[0], state[1], state[2]];
    let v = [state[3], state[4], state[5]];
    // Resolve ΔV and execution covariance into ECI.
    let (dv_eci, q_eci) = match man.frame {
        ManeuverFrame::Eci => (man.dv, man.exec_cov),
        ManeuverFrame::Lvlh => {
            let m = lvlh_to_eci(r, v);
            // Q_eci = M · exec_cov · Mᵀ.
            let mut q = [[0.0f64; 3]; 3];
            for (i, qi) in q.iter_mut().enumerate() {
                for (j, qij) in qi.iter_mut().enumerate() {
                    let mut s = 0.0;
                    for (a, ma) in man.exec_cov.iter().enumerate() {
                        for (b, &mab) in ma.iter().enumerate() {
                            // M[i][a] · exec[a][b] · M[j][b], with M[row=i] = (col a)[i].
                            s += m[a][i] * mab * m[b][j];
                        }
                    }
                    *qij = s;
                }
            }
            (rotate(&m, man.dv), q)
        }
    };
    let mut out = state;
    out[3] += dv_eci[0];
    out[4] += dv_eci[1];
    out[5] += dv_eci[2];
    let mut cov_after = cov;
    for i in 0..3 {
        for j in 0..3 {
            cov_after[3 + i][3 + j] += q_eci[i][j];
        }
    }
    (out, cov_after)
}

// ===========================================================================
// 2. Finite-burn maneuver vs the Tsiolkovsky rocket equation.
// ===========================================================================

/// A constant-thrust burn over a fixed arc, integrated with mass as a state variable.
#[derive(Clone, Copy, Debug)]
pub struct FiniteBurn {
    /// Engine thrust (N).
    pub thrust_n: f64,
    /// Specific impulse (s); effective exhaust velocity is `isp_s · G0`.
    pub isp_s: f64,
    /// Initial wet mass (kg).
    pub m0_kg: f64,
    /// Burn duration (s).
    pub burn_s: f64,
    /// Thrust direction (need not be normalized; only its direction is used).
    pub dir: Vec3,
}

/// Outcome of integrating a [`FiniteBurn`].
#[derive(Clone, Copy, Debug, Serialize)]
pub struct FiniteBurnResult {
    /// Achieved |ΔV| (m/s) from the field-free numerical integration.
    pub dv_ms: f64,
    /// Final (dry-of-burn) mass (kg).
    pub mf_kg: f64,
    /// Closed-form Tsiolkovsky ΔV = `Isp·g₀·ln(m0/mf)` (m/s).
    pub tsiolkovsky_ms: f64,
    /// Relative error |dv − tsiolkovsky| / tsiolkovsky.
    pub rel_err: f64,
}

/// Tsiolkovsky rocket equation: ΔV = `Isp · g₀ · ln(m0 / mf)` (m/s).
pub fn tsiolkovsky(isp_s: f64, m0_kg: f64, mf_kg: f64) -> f64 {
    isp_s * G0 * (m0_kg / mf_kg).pln()
}

/// Integrate a finite burn in a field-free setting (isolating the rocket-equation check from
/// gravity losses) with `steps` RK4 sub-steps. The coupled ODE on `[v(3), m]` is
/// `v̇ = (T/m)·d̂`, `ṁ = −T/(Isp·g₀)`; the achieved |ΔV| is compared to Tsiolkovsky.
///
/// Returns an error if the burn would exhaust the mass (`mf ≤ 0`).
pub fn integrate_finite_burn(burn: &FiniteBurn, steps: usize) -> Result<FiniteBurnResult, String> {
    let mdot = burn.thrust_n / (burn.isp_s * G0);
    let mf = burn.m0_kg - mdot * burn.burn_s;
    if mf <= 0.0 {
        return Err(format!(
            "burn exhausts mass: m0={} kg, ṁ={} kg/s, burn={} s ⇒ mf={} kg",
            burn.m0_kg, mdot, burn.burn_s, mf
        ));
    }
    let d = unit(burn.dir);
    let h = burn.burn_s / steps as f64;
    // y = [vx, vy, vz, m]; field-free, so v̇ = (T/m)·d̂, ṁ = −mdot.
    let deriv = |_t: f64, y: &[f64]| -> Vec<f64> {
        let m = y[3];
        let a = burn.thrust_n / m;
        vec![a * d[0], a * d[1], a * d[2], -mdot]
    };
    let mut y = vec![0.0, 0.0, 0.0, burn.m0_kg];
    let mut t = 0.0;
    for _ in 0..steps {
        y = crate::integrator::rk4_step(&deriv, t, &y, h);
        t += h;
    }
    let dv_ms = norm([y[0], y[1], y[2]]);
    let tsiol = tsiolkovsky(burn.isp_s, burn.m0_kg, mf);
    let rel_err = (dv_ms - tsiol).abs() / tsiol;
    Ok(FiniteBurnResult {
        dv_ms,
        mf_kg: mf,
        tsiolkovsky_ms: tsiol,
        rel_err,
    })
}

// ===========================================================================
// 3. Exact universal-variable Kepler propagator (two-body truth for validation).
// ===========================================================================

/// Stumpff functions `c2(ψ)`, `c3(ψ)` (series near zero, closed form otherwise).
fn stumpff(psi: f64) -> (f64, f64) {
    if psi > 1e-6 {
        let s = psi.sqrt();
        ((1.0 - s.pcos()) / psi, (s - s.psin()) / (psi * s))
    } else if psi < -1e-6 {
        let s = (-psi).sqrt();
        ((s.pcosh() - 1.0) / (-psi), (s.psinh() - s) / ((-psi) * s))
    } else {
        // Series: c2 = 1/2 − ψ/24 + ψ²/720…, c3 = 1/6 − ψ/120 + ψ²/5040…
        (
            0.5 - psi / 24.0 + psi * psi / 720.0,
            1.0 / 6.0 - psi / 120.0 + psi * psi / 5040.0,
        )
    }
}

/// Propagate a two-body state `(r0, v0)` forward by `dt` under gravity parameter `mu`
/// using the universal-variable (Stumpff f/g) formulation. Exact two-body to convergence —
/// the reference truth the Lambert outputs are round-tripped against.
pub fn kepler_universal(r0: Vec3, v0: Vec3, dt: f64, mu: f64) -> (Vec3, Vec3) {
    let sqrt_mu = mu.sqrt();
    let r0n = norm(r0);
    let v0n = norm(v0);
    let rv0 = dot(r0, v0); // r0·v0 (NOT divided by |r0| — the universal-variable term needs r0·v0/√μ)
    let alpha = 2.0 / r0n - v0n * v0n / mu; // 1/a
                                            // Initial χ guess.
    let mut chi = if alpha > 1e-9 {
        sqrt_mu * dt * alpha
    } else if alpha < -1e-9 {
        let a = 1.0 / alpha;
        dt.signum()
            * (-a).sqrt()
            * ((-2.0 * mu * alpha * dt)
                / (dot(r0, v0) + dt.signum() * (-mu * a).sqrt() * (1.0 - r0n * alpha)))
                .pln()
    } else {
        sqrt_mu * dt / r0n
    };
    for _ in 0..100 {
        let psi = chi * chi * alpha;
        let (c2, c3) = stumpff(psi);
        let r = chi * chi * c2 + (rv0 / sqrt_mu) * chi * (1.0 - psi * c3) + r0n * (1.0 - psi * c2);
        let f = sqrt_mu * dt
            - chi * chi * chi * c3
            - (rv0 / sqrt_mu) * chi * chi * c2
            - r0n * chi * (1.0 - psi * c3);
        let dchi = f / r;
        chi += dchi;
        if dchi.abs() < 1e-10 {
            break;
        }
    }
    let psi = chi * chi * alpha;
    let (c2, c3) = stumpff(psi);
    let f = 1.0 - chi * chi / r0n * c2;
    let g = dt - chi * chi * chi / sqrt_mu * c3;
    let r_vec = add(scale(r0, f), scale(v0, g));
    let rn = norm(r_vec);
    let fdot = sqrt_mu / (rn * r0n) * chi * (psi * c3 - 1.0);
    let gdot = 1.0 - chi * chi / rn * c2;
    let v_vec = add(scale(r0, fdot), scale(v0, gdot));
    (r_vec, v_vec)
}

// ===========================================================================
// 4. Izzo-2015 single-revolution Lambert solver.
// ===========================================================================

/// Gauss hypergeometric `₂F₁(3, 1; 5/2; x)` by its power series (Izzo's `hyp2f1b`),
/// used for the time-of-flight series near `x = 1`.
fn hyp2f1b(x: f64) -> f64 {
    if x >= 1.0 {
        return f64::INFINITY;
    }
    let mut res = 1.0_f64;
    let mut term = 1.0_f64;
    let mut i = 0.0_f64;
    loop {
        term *= (3.0 + i) * (1.0 + i) / (2.5 + i) * x / (i + 1.0);
        let newres = res + term;
        if newres == res {
            return newres;
        }
        res = newres;
        i += 1.0;
    }
}

fn compute_y(x: f64, ll: f64) -> f64 {
    (1.0 - ll * ll * (1.0 - x * x)).sqrt()
}

fn compute_psi(x: f64, y: f64, ll: f64) -> f64 {
    if (-1.0..1.0).contains(&x) {
        // Elliptic.
        (x * y + ll * (1.0 - x * x)).pacos()
    } else if x > 1.0 {
        // Hyperbolic.
        ((y - x * ll) * (x * x - 1.0).sqrt()).pasinh()
    } else {
        0.0
    }
}

/// `tof(x) − T0` for the (non-dimensional) Lambert time equation, single revolution.
fn tof_equation_y(x: f64, y: f64, t0: f64, ll: f64) -> f64 {
    let tof = if (0.6_f64).sqrt() < x && x < (1.4_f64).sqrt() {
        let eta = y - ll * x;
        let s1 = 0.5 * (1.0 - ll - x * eta);
        let q = 4.0 / 3.0 * hyp2f1b(s1);
        0.5 * (eta * eta * eta * q + 4.0 * ll * eta)
    } else {
        let psi = compute_psi(x, y, ll);
        (psi / (1.0 - x * x).abs().sqrt() - x + ll * y) / (1.0 - x * x)
    };
    tof - t0
}

fn tof_eq_p(x: f64, y: f64, tof: f64, ll: f64) -> f64 {
    (3.0 * tof * x - 2.0 + 2.0 * ll * ll * ll * x / y) / (1.0 - x * x)
}
fn tof_eq_p2(x: f64, y: f64, tof: f64, dt: f64, ll: f64) -> f64 {
    (3.0 * tof + 5.0 * x * dt + 2.0 * (1.0 - ll * ll) * ll * ll * ll / (y * y * y)) / (1.0 - x * x)
}
fn tof_eq_p3(x: f64, y: f64, _tof: f64, dt: f64, ddt: f64, ll: f64) -> f64 {
    (7.0 * x * ddt + 8.0 * dt - 6.0 * (1.0 - ll * ll) * ll.ppowi(5) * x / y.ppowi(5))
        / (1.0 - x * x)
}

/// Householder (cubic-order) iteration for the single-rev root `x` of `tof(x) = T`.
fn householder(x0: f64, big_t: f64, ll: f64) -> f64 {
    let mut x = x0;
    for _ in 0..35 {
        let y = compute_y(x, ll);
        let fval = tof_equation_y(x, y, big_t, ll);
        let tof = fval + big_t;
        let fp = tof_eq_p(x, y, tof, ll);
        let fpp = tof_eq_p2(x, y, tof, fp, ll);
        let fppp = tof_eq_p3(x, y, tof, fp, fpp, ll);
        let denom = fp * (fp * fp - fval * fpp) + fppp * fval * fval / 6.0;
        let dx = fval * (fp * fp - fval * fpp / 2.0) / denom;
        x -= dx;
        if dx.abs() < 1e-12 {
            break;
        }
    }
    x
}

/// Solve Lambert's problem (single revolution, `M = 0`) by Izzo's 2015 method.
///
/// Given two position vectors `r1`, `r2`, the time of flight `tof` (s), gravity parameter `mu`,
/// and a `prograde` flag (transfer with `+z` angular momentum), returns the departure and
/// arrival velocities `(v1, v2)` of the connecting two-body arc. Errors on degenerate geometry
/// (collinear `r1`, `r2` ⇒ undefined transfer plane) or non-positive time of flight.
pub fn lambert(
    r1: Vec3,
    r2: Vec3,
    tof: f64,
    mu: f64,
    prograde: bool,
) -> Result<(Vec3, Vec3), String> {
    if tof <= 0.0 {
        return Err("time of flight must be positive".into());
    }
    let r1n = norm(r1);
    let r2n = norm(r2);
    let c_vec = sub(r2, r1);
    let c = norm(c_vec);
    let s = 0.5 * (r1n + r2n + c);
    let ir1 = unit(r1);
    let ir2 = unit(r2);
    let cross12 = cross(ir1, ir2);
    let cross_n = norm(cross12);
    if cross_n < 1e-12 {
        return Err("collinear endpoints: transfer plane is undefined (≈0°/180°)".into());
    }
    let ih = scale(cross12, 1.0 / cross_n);
    let mut ll = (1.0 - (c / s).min(1.0)).sqrt();
    let (mut it1, mut it2);
    if ih[2] < 0.0 {
        ll = -ll;
        it1 = cross(ir1, ih);
        it2 = cross(ir2, ih);
    } else {
        it1 = cross(ih, ir1);
        it2 = cross(ih, ir2);
    }
    if !prograde {
        ll = -ll;
        it1 = scale(it1, -1.0);
        it2 = scale(it2, -1.0);
    }
    it1 = unit(it1);
    it2 = unit(it2);

    let big_t = (2.0 * mu / (s * s * s)).sqrt() * tof;

    // Single-revolution initial guess (Izzo `_initial_guess`, M = 0 branch).
    let t0 = ll.pacos() + ll * (1.0 - ll * ll).sqrt(); // T(x = 0)
    let t1 = 2.0 / 3.0 * (1.0 - ll * ll * ll); // T(x = 1), parabola
    let x0 = if big_t >= t0 {
        (t0 / big_t).ppowf(2.0 / 3.0) - 1.0
    } else if big_t < t1 {
        2.5 * t1 / big_t * (t1 - big_t) / (1.0 - ll.ppowi(5)) + 1.0
    } else {
        (t0 / big_t).ppowf((t1 / t0).plog2()) - 1.0
    };
    let x = householder(x0, big_t, ll);
    let y = compute_y(x, ll);

    let gamma = (mu * s / 2.0).sqrt();
    let rho = (r1n - r2n) / c;
    let sigma = (1.0 - rho * rho).sqrt();
    let vr1 = gamma * ((ll * y - x) - rho * (ll * y + x)) / r1n;
    let vr2 = -gamma * ((ll * y - x) + rho * (ll * y + x)) / r2n;
    let vt1 = gamma * sigma * (y + ll * x) / r1n;
    let vt2 = gamma * sigma * (y + ll * x) / r2n;
    let v1 = add(scale(ir1, vr1), scale(it1, vt1));
    let v2 = add(scale(ir2, vr2), scale(it2, vt2));
    Ok((v1, v2))
}

// ===========================================================================
// 4b. Izzo-2015 multi-revolution Lambert branches (M >= 1).
// ===========================================================================

/// Which of the two `M`-revolution Lambert solutions (Izzo 2015, §3): the **left** branch starts
/// its Householder iteration from Izzo's left initial guess `x0l`, the **right** branch from
/// `x0r`. The names and the order follow Izzo's paper and his reference implementation.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Serialize)]
#[serde(rename_all = "lowercase")]
pub enum LambertBranch {
    /// The solution reached from the left initial guess.
    Left,
    /// The solution reached from the right initial guess.
    Right,
}

/// One solution of Lambert's problem.
#[derive(Clone, Copy, Debug, PartialEq, Serialize)]
pub struct LambertSolution {
    /// Complete revolutions `M` before arrival.
    pub revs: u32,
    /// `None` for `M = 0` (unique); the branch for `M >= 1`.
    pub branch: Option<LambertBranch>,
    /// Departure velocity (same units as `mu` and the positions imply; m/s in SI).
    pub v1: Vec3,
    /// Arrival velocity.
    pub v2: Vec3,
}

/// The non-dimensional Lambert geometry of Izzo 2015, §2: `λ`, the tangential unit vectors and
/// the scale factors. Same construction as [`lambert`].
struct LambertGeom {
    ll: f64,
    ir1: Vec3,
    ir2: Vec3,
    it1: Vec3,
    it2: Vec3,
    r1n: f64,
    r2n: f64,
    c: f64,
    s: f64,
}

fn lambert_geom(r1: Vec3, r2: Vec3, prograde: bool) -> Result<LambertGeom, String> {
    let r1n = norm(r1);
    let r2n = norm(r2);
    let c = norm(sub(r2, r1));
    let s = 0.5 * (r1n + r2n + c);
    let ir1 = unit(r1);
    let ir2 = unit(r2);
    let cross12 = cross(ir1, ir2);
    let cross_n = norm(cross12);
    if cross_n < 1e-12 {
        return Err("collinear endpoints: transfer plane is undefined (≈0°/180°)".into());
    }
    let ih = scale(cross12, 1.0 / cross_n);
    let mut ll = (1.0 - (c / s).min(1.0)).sqrt();
    let (mut it1, mut it2);
    if ih[2] < 0.0 {
        ll = -ll;
        it1 = cross(ir1, ih);
        it2 = cross(ir2, ih);
    } else {
        it1 = cross(ih, ir1);
        it2 = cross(ih, ir2);
    }
    if !prograde {
        ll = -ll;
        it1 = scale(it1, -1.0);
        it2 = scale(it2, -1.0);
    }
    Ok(LambertGeom {
        ll,
        ir1,
        ir2,
        it1: unit(it1),
        it2: unit(it2),
        r1n,
        r2n,
        c,
        s,
    })
}

/// Non-dimensional time of flight `T(x)` on the elliptic `M`-revolution family (`|x| < 1`):
/// `T = [(ψ + Mπ)/√(1 − x²) − x + λy] / (1 − x²)` (Izzo 2015, eq. 18 with the `Mπ` term).
fn tof_multirev(x: f64, ll: f64, m: u32) -> f64 {
    let y = compute_y(x, ll);
    let psi = compute_psi(x, y, ll);
    let om = 1.0 - x * x;
    ((psi + m as f64 * PI) / om.sqrt() - x + ll * y) / om
}

/// Householder (cubic-order) iteration on `T(x) = big_t` along the `M`-revolution family.
fn householder_multirev(x0: f64, big_t: f64, ll: f64, m: u32) -> f64 {
    let mut x = x0;
    for _ in 0..60 {
        let y = compute_y(x, ll);
        let tof = tof_multirev(x, ll, m);
        let fval = tof - big_t;
        let fp = tof_eq_p(x, y, tof, ll);
        let fpp = tof_eq_p2(x, y, tof, fp, ll);
        let fppp = tof_eq_p3(x, y, tof, fp, fpp, ll);
        let denom = fp * (fp * fp - fval * fpp) + fppp * fval * fval / 6.0;
        let dx = fval * (fp * fp - fval * fpp / 2.0) / denom;
        // Keep the iterate inside the ellipse family.
        let next = (x - dx).clamp(-1.0 + 1e-15, 1.0 - 1e-15);
        let step = (next - x).abs();
        x = next;
        if step < 1e-14 {
            break;
        }
    }
    x
}

/// The minimum non-dimensional time of flight `T_min(M)` of the `M`-revolution family and its
/// `x`: Halley iteration on `dT/dx = 0` from `x = 0.1` (Izzo 2015, §3).
fn tmin_multirev(ll: f64, m: u32) -> (f64, f64) {
    if ll.abs() == 1.0 {
        return (0.0, tof_multirev(0.0, ll, m));
    }
    let mut x = 0.1_f64;
    for _ in 0..60 {
        let y = compute_y(x, ll);
        let tof = tof_multirev(x, ll, m);
        let d1 = tof_eq_p(x, y, tof, ll);
        let d2 = tof_eq_p2(x, y, tof, d1, ll);
        let d3 = tof_eq_p3(x, y, tof, d1, d2, ll);
        let den = 2.0 * d2 * d2 - d1 * d3;
        if den == 0.0 {
            break;
        }
        let next = (x - 2.0 * d1 * d2 / den).clamp(-1.0 + 1e-15, 1.0 - 1e-15);
        let step = (next - x).abs();
        x = next;
        if step < 1e-14 {
            break;
        }
    }
    (x, tof_multirev(x, ll, m))
}

/// Terminal velocities from a converged `x` (Izzo 2015, §2, Algorithm 1 last step).
fn lambert_velocities(g: &LambertGeom, x: f64, mu: f64) -> (Vec3, Vec3) {
    let ll = g.ll;
    let y = compute_y(x, ll);
    let gamma = (mu * g.s / 2.0).sqrt();
    let rho = (g.r1n - g.r2n) / g.c;
    let sigma = (1.0 - rho * rho).sqrt();
    let vr1 = gamma * ((ll * y - x) - rho * (ll * y + x)) / g.r1n;
    let vr2 = -gamma * ((ll * y - x) + rho * (ll * y + x)) / g.r2n;
    let vt1 = gamma * sigma * (y + ll * x) / g.r1n;
    let vt2 = gamma * sigma * (y + ll * x) / g.r2n;
    (
        add(scale(g.ir1, vr1), scale(g.it1, vt1)),
        add(scale(g.ir2, vr2), scale(g.it2, vt2)),
    )
}

/// The largest number of complete revolutions `M <= max_revs` for which a Lambert solution
/// exists: `floor(T/π)` cropped to `max_revs`, lowered by one when `T` is below that family's
/// minimum time of flight (Izzo 2015, §3).
pub fn lambert_max_revs(
    r1: Vec3,
    r2: Vec3,
    tof: f64,
    mu: f64,
    prograde: bool,
    max_revs: u32,
) -> Result<u32, String> {
    if tof <= 0.0 {
        return Err("time of flight must be positive".into());
    }
    let g = lambert_geom(r1, r2, prograde)?;
    let big_t = (2.0 * mu / (g.s * g.s * g.s)).sqrt() * tof;
    Ok(max_revs_of(&g, big_t, max_revs))
}

fn max_revs_of(g: &LambertGeom, big_t: f64, max_revs: u32) -> u32 {
    let ll = g.ll;
    let mut n = ((big_t / PI).floor() as u32).min(max_revs);
    let t00 = ll.pacos() + ll * (1.0 - ll * ll).sqrt();
    if n > 0 && big_t < t00 + n as f64 * PI {
        let (_, t_min) = tmin_multirev(ll, n);
        if t_min > big_t {
            n -= 1;
        }
    }
    n
}

/// Every solution of Lambert's problem with at most `max_revs` complete revolutions, by Izzo's
/// 2015 method: the single-revolution solution first (exactly [`lambert`]), then for each
/// `M = 1 ..= Nmax` the left and the right branch. `Nmax` is [`lambert_max_revs`]; with
/// `max_revs = 0` this is `[lambert(..)]`.
///
/// Near a family's minimum time of flight the two branches merge and the root is
/// ill-conditioned; there the velocities lose digits (see the oracle test's bars).
pub fn lambert_multirev(
    r1: Vec3,
    r2: Vec3,
    tof: f64,
    mu: f64,
    prograde: bool,
    max_revs: u32,
) -> Result<Vec<LambertSolution>, String> {
    let (v1, v2) = lambert(r1, r2, tof, mu, prograde)?;
    let mut out = vec![LambertSolution {
        revs: 0,
        branch: None,
        v1,
        v2,
    }];
    let g = lambert_geom(r1, r2, prograde)?;
    let big_t = (2.0 * mu / (g.s * g.s * g.s)).sqrt() * tof;
    let n = max_revs_of(&g, big_t, max_revs);
    for m in 1..=n {
        let mf = m as f64;
        let kl = ((mf * PI + PI) / (8.0 * big_t)).ppowf(2.0 / 3.0);
        let kr = ((8.0 * big_t) / (mf * PI)).ppowf(2.0 / 3.0);
        for (branch, k) in [(LambertBranch::Left, kl), (LambertBranch::Right, kr)] {
            let x = householder_multirev((k - 1.0) / (k + 1.0), big_t, g.ll, m);
            let (v1, v2) = lambert_velocities(&g, x, mu);
            out.push(LambertSolution {
                revs: m,
                branch: Some(branch),
                v1,
                v2,
            });
        }
    }
    Ok(out)
}

/// One `M`-revolution Lambert solution on the given branch (`M >= 1`), or an error when no
/// `M`-revolution transfer exists for this time of flight. For `M = 0` use [`lambert`].
pub fn lambert_revs(
    r1: Vec3,
    r2: Vec3,
    tof: f64,
    mu: f64,
    prograde: bool,
    revs: u32,
    branch: LambertBranch,
) -> Result<(Vec3, Vec3), String> {
    if revs == 0 {
        return Err("lambert_revs is for M >= 1; use lambert for M = 0".into());
    }
    let all = lambert_multirev(r1, r2, tof, mu, prograde, revs)?;
    all.iter()
        .find(|s| s.revs == revs && s.branch == Some(branch))
        .map(|s| (s.v1, s.v2))
        .ok_or_else(|| {
            format!("no {revs}-revolution transfer: the time of flight is below its minimum")
        })
}

// ===========================================================================
// 5. Porkchop sweep over a synthetic coplanar-circular heliocentric system.
// ===========================================================================

/// A body on a coplanar circular heliocentric orbit (synthetic, so the transfer optimum is
/// analytic). `phase0_rad` is its true longitude at epoch `t = 0`.
#[derive(Clone, Copy, Debug)]
pub struct CircularCoplanarBody {
    /// Orbit radius (m).
    pub radius_m: f64,
    /// True longitude at `t = 0` (rad).
    pub phase0_rad: f64,
    /// Central gravity parameter (m³/s²).
    pub mu_central: f64,
}

impl CircularCoplanarBody {
    /// Mean angular rate `n = √(μ/r³)` (rad/s).
    pub fn angular_rate(&self) -> f64 {
        (self.mu_central / self.radius_m.ppowi(3)).sqrt()
    }
    /// Heliocentric position (m) at time `t` (s).
    pub fn position(&self, t: f64) -> Vec3 {
        let th = self.phase0_rad + self.angular_rate() * t;
        [self.radius_m * th.pcos(), self.radius_m * th.psin(), 0.0]
    }
    /// Heliocentric velocity (m/s) at time `t` (s).
    pub fn velocity(&self, t: f64) -> Vec3 {
        let n = self.angular_rate();
        let th = self.phase0_rad + n * t;
        let speed = self.radius_m * n;
        [-speed * th.psin(), speed * th.pcos(), 0.0]
    }
}

/// A porkchop grid: departure C3 (km²/s²) and arrival V∞ (km/s) over launch × arrival epochs.
/// Degenerate cells (non-positive TOF or near-collinear ≈0°/180° geometry) are stored as NaN.
#[derive(Clone, Debug, Serialize)]
pub struct PorkchopGrid {
    /// Departure epochs (s).
    pub dep_epochs_s: Vec<f64>,
    /// Arrival epochs (s).
    pub arr_epochs_s: Vec<f64>,
    /// `c3_km2s2[i][j]` = departure C3 for `dep_epochs_s[i]` → `arr_epochs_s[j]`.
    pub c3_km2s2: Vec<Vec<f64>>,
    /// `vinf_arr_kms[i][j]` = arrival hyperbolic-excess speed for the same cell.
    pub vinf_arr_kms: Vec<Vec<f64>>,
}

impl PorkchopGrid {
    /// Serialize the grid (epoch axes + the two 2-D arrays) for browser contour rendering.
    pub fn to_json(&self) -> String {
        // `PorkchopGrid` is only `Vec<f64>` / `Vec<Vec<f64>>` fields — no maps and no
        // fallible custom `Serialize` impls — so JSON serialisation cannot fail. (A
        // non-finite C3/v∞ entry serialises to `null`; it does not error.)
        serde_json::to_string_pretty(self)
            .expect("PorkchopGrid (only Vec<f64>/Vec<Vec<f64>> fields) always serialises")
    }
    /// The minimum finite C3 over the grid and its `(dep_index, arr_index)`.
    pub fn min_c3(&self) -> Option<(f64, usize, usize)> {
        let mut best: Option<(f64, usize, usize)> = None;
        for (i, row) in self.c3_km2s2.iter().enumerate() {
            for (j, &c3) in row.iter().enumerate() {
                let better = match best {
                    None => true,
                    Some((b, _, _)) => c3 < b,
                };
                if c3.is_finite() && better {
                    best = Some((c3, i, j));
                }
            }
        }
        best
    }
}

/// Sweep a launch-epoch × arrival-epoch grid, solving Lambert for each cell and recording the
/// departure C3 (km²/s²) and arrival V∞ (km/s). Cells within `0.5°` of a 0°/180° transfer (or
/// with TOF ≤ 0) are marked NaN — the transfer plane is ill-conditioned there.
pub fn porkchop(
    dep: &CircularCoplanarBody,
    arr: &CircularCoplanarBody,
    dep_epochs_s: &[f64],
    arr_epochs_s: &[f64],
    mu_helio: f64,
) -> PorkchopGrid {
    let collinear_guard = (0.5_f64).to_radians();
    let mut c3 = vec![vec![f64::NAN; arr_epochs_s.len()]; dep_epochs_s.len()];
    let mut vinf = vec![vec![f64::NAN; arr_epochs_s.len()]; dep_epochs_s.len()];
    for (i, &td) in dep_epochs_s.iter().enumerate() {
        for (j, &ta) in arr_epochs_s.iter().enumerate() {
            let tof = ta - td;
            if tof <= 0.0 {
                continue;
            }
            let r1 = dep.position(td);
            let r2 = arr.position(ta);
            // Skip near-collinear geometry.
            let cosang = (dot(r1, r2) / (norm(r1) * norm(r2))).clamp(-1.0, 1.0);
            let ang = cosang.pacos();
            if ang < collinear_guard || (PI - ang) < collinear_guard {
                continue;
            }
            if let Ok((v1, v2)) = lambert(r1, r2, tof, mu_helio, true) {
                let vinf_dep = sub(v1, dep.velocity(td));
                let vinf_arr = sub(v2, arr.velocity(ta));
                c3[i][j] = dot(vinf_dep, vinf_dep) / 1.0e6; // m²/s² → km²/s²
                vinf[i][j] = norm(vinf_arr) / 1000.0; // m/s → km/s
            }
        }
    }
    PorkchopGrid {
        dep_epochs_s: dep_epochs_s.to_vec(),
        arr_epochs_s: arr_epochs_s.to_vec(),
        c3_km2s2: c3,
        vinf_arr_kms: vinf,
    }
}

/// The porkchop sweep of [`porkchop`] on a JPL planetary ephemeris: the departure and arrival
/// states are the heliocentric J2000 states of the NAIF bodies `dep_body` and `arr_body` (for
/// example 399 for the Earth and 4 for the Mars barycentre) read from `kernel` (a DE440 SPK) at
/// the TDB Julian dates `dep_jd_tdb` and `arr_jd_tdb`, with heliocentric parameter `mu_sun`.
/// Same single-revolution prograde Lambert, the same 0.5° collinearity guard, C3 in km²/s² and
/// arrival V∞ in km/s; the epoch axes of the grid are seconds from the first departure date.
/// Checked against SPICE states and pykep's Lambert in `tests/porkchop_de440_oracle.rs`.
///
/// # Errors
/// When the kernel does not cover an epoch or lacks a body.
pub fn porkchop_kernel(
    kernel: &crate::ephem_provider::KernelEphemeris,
    dep_body: i32,
    arr_body: i32,
    dep_jd_tdb: &[f64],
    arr_jd_tdb: &[f64],
    mu_sun: f64,
) -> Result<PorkchopGrid, String> {
    use crate::ephem_provider::naif_code::SUN;
    let st = |b: i32, jd: f64| kernel.state(b, SUN, crate::jd2::Jd2::new(jd));
    let dep: Vec<_> = dep_jd_tdb
        .iter()
        .map(|&j| st(dep_body, j))
        .collect::<Result<_, _>>()?;
    let arr: Vec<_> = arr_jd_tdb
        .iter()
        .map(|&j| st(arr_body, j))
        .collect::<Result<_, _>>()?;
    let t0 = dep_jd_tdb.first().copied().unwrap_or(0.0);
    let collinear_guard = (0.5_f64).to_radians();
    let mut c3 = vec![vec![f64::NAN; arr_jd_tdb.len()]; dep_jd_tdb.len()];
    let mut vinf = vec![vec![f64::NAN; arr_jd_tdb.len()]; dep_jd_tdb.len()];
    for (i, &jd) in dep_jd_tdb.iter().enumerate() {
        for (j, &ja) in arr_jd_tdb.iter().enumerate() {
            let tof = (ja - jd) * 86_400.0;
            if tof <= 0.0 {
                continue;
            }
            let (r1, r2) = (dep[i][0], arr[j][0]);
            let cosang = (dot(r1, r2) / (norm(r1) * norm(r2))).clamp(-1.0, 1.0);
            let ang = cosang.pacos();
            if ang < collinear_guard || (PI - ang) < collinear_guard {
                continue;
            }
            if let Ok((v1, v2)) = lambert(r1, r2, tof, mu_sun, true) {
                let vd = sub(v1, dep[i][1]);
                let va = sub(v2, arr[j][1]);
                c3[i][j] = dot(vd, vd) / 1.0e6;
                vinf[i][j] = norm(va) / 1000.0;
            }
        }
    }
    Ok(PorkchopGrid {
        dep_epochs_s: dep_jd_tdb.iter().map(|j| (j - t0) * 86_400.0).collect(),
        arr_epochs_s: arr_jd_tdb.iter().map(|j| (j - t0) * 86_400.0).collect(),
        c3_km2s2: c3,
        vinf_arr_kms: vinf,
    })
}

/// What [`optimise_transfer_kernel`] minimises.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum TransferObjective {
    /// Departure C3 `|v1 − v_dep|²` (km²/s²).
    MinC3,
    /// Arrival hyperbolic excess speed `|v2 − v_arr|` (km/s).
    MinArrivalVinf,
}

/// Transfer type: the heliocentric transfer angle below (Type I) or above (Type II) 180°,
/// measured in the prograde sense about the frame's +Z axis (the convention of the prograde
/// Lambert solver), so Type I means `(r1 × r2)_z > 0`.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum TransferType {
    /// Transfer angle below 180°.
    I,
    /// Transfer angle above 180°.
    II,
}

/// The optimum found by [`optimise_transfer_kernel`].
#[derive(Clone, Copy, Debug)]
pub struct TransferOptimum {
    /// Departure epoch (JD, TDB).
    pub dep_jd_tdb: f64,
    /// Arrival epoch (JD, TDB).
    pub arr_jd_tdb: f64,
    /// Departure C3 (km²/s²).
    pub c3_km2s2: f64,
    /// Arrival hyperbolic excess speed (km/s).
    pub vinf_arr_kms: f64,
    /// Prograde transfer angle (deg, 0..360).
    pub transfer_angle_deg: f64,
    /// Nelder-Mead iterations used.
    pub iterations: usize,
}

/// One ballistic two-impulse transfer between the states of two NAIF bodies on the kernel:
/// `(C3 km²/s², arrival V∞ km/s, prograde transfer angle deg)`, or `None` for a non-positive
/// time of flight, a geometry within the 0.5° collinearity guard, the other transfer type, or a
/// Lambert failure.
#[allow(clippy::too_many_arguments)]
fn transfer_eval(
    kernel: &crate::ephem_provider::KernelEphemeris,
    dep_body: i32,
    arr_body: i32,
    dep_jd: f64,
    arr_jd: f64,
    ttype: TransferType,
    mu_sun: f64,
) -> Option<(f64, f64, f64)> {
    use crate::ephem_provider::naif_code::SUN;
    let tof = (arr_jd - dep_jd) * 86_400.0;
    if tof.is_nan() || tof <= 0.0 {
        return None;
    }
    let d = kernel
        .state(dep_body, SUN, crate::jd2::Jd2::new(dep_jd))
        .ok()?;
    let a = kernel
        .state(arr_body, SUN, crate::jd2::Jd2::new(arr_jd))
        .ok()?;
    let (r1, r2) = (d[0], a[0]);
    let cosang = (dot(r1, r2) / (norm(r1) * norm(r2))).clamp(-1.0, 1.0);
    let ang = cosang.pacos();
    let guard = (0.5_f64).to_radians();
    if ang < guard || (PI - ang) < guard {
        return None;
    }
    let short = cross(r1, r2)[2] > 0.0;
    if short != (ttype == TransferType::I) {
        return None;
    }
    let theta = if short { ang } else { 2.0 * PI - ang };
    let (v1, v2) = lambert(r1, r2, tof, mu_sun, true).ok()?;
    let vd = sub(v1, d[1]);
    let va = sub(v2, a[1]);
    Some((dot(vd, vd) / 1.0e6, norm(va) / 1000.0, theta.to_degrees()))
}

/// The departure and arrival dates of a single-revolution ballistic two-impulse transfer that
/// minimise `objective` within one transfer type, as a continuous optimum (CAP-20c): a grid
/// over both windows at `grid_step_days`, then 2-D Nelder-Mead on (departure, arrival) in days
/// from the best cell (initial steps 1 day; stop when the simplex spans < 1e-7 day and the
/// objective spread is < 1e-13 relative, or after 4000 iterations). Candidates of the other
/// type, with a non-positive time of flight, within the collinearity guard or where Lambert
/// fails score +infinity. States come from the DE440 SPK as in [`porkchop_kernel`].
///
/// A date optimiser for one ballistic leg, not a mission design: no launch-vehicle, parking
/// orbit, deep-space manoeuvre or multi-revolution options.
#[allow(clippy::too_many_arguments)]
pub fn optimise_transfer_kernel(
    kernel: &crate::ephem_provider::KernelEphemeris,
    dep_body: i32,
    arr_body: i32,
    dep_window_jd_tdb: (f64, f64),
    arr_window_jd_tdb: (f64, f64),
    objective: TransferObjective,
    ttype: TransferType,
    mu_sun: f64,
    grid_step_days: f64,
) -> Result<TransferOptimum, String> {
    if !(grid_step_days.is_finite() && grid_step_days > 0.0) {
        return Err(format!(
            "grid_step_days must be positive, got {grid_step_days}"
        ));
    }
    let (d0, d1) = dep_window_jd_tdb;
    let (a0, a1) = arr_window_jd_tdb;
    if !(d1 >= d0 && a1 >= a0) {
        return Err("each window must run forward (start <= end)".into());
    }
    let pick = |e: (f64, f64, f64)| match objective {
        TransferObjective::MinC3 => e.0,
        TransferObjective::MinArrivalVinf => e.1,
    };
    let f = |x: [f64; 2]| -> f64 {
        transfer_eval(kernel, dep_body, arr_body, x[0], x[1], ttype, mu_sun)
            .map(pick)
            .unwrap_or(f64::INFINITY)
    };
    // Grid.
    let steps = |lo: f64, hi: f64| ((hi - lo) / grid_step_days + 1e-9).floor() as usize + 1;
    let mut best = ([f64::NAN; 2], f64::INFINITY);
    for i in 0..steps(d0, d1) {
        let dj = d0 + i as f64 * grid_step_days;
        for j in 0..steps(a0, a1) {
            let aj = a0 + j as f64 * grid_step_days;
            let v = f([dj, aj]);
            if v < best.1 {
                best = ([dj, aj], v);
            }
        }
    }
    if !best.1.is_finite() {
        return Err("no transfer of the requested type in the windows".into());
    }
    // Nelder-Mead (reflection 1, expansion 2, contraction 0.5, shrink 0.5).
    let x0 = best.0;
    let mut simplex = [
        (x0, best.1),
        ([x0[0] + 1.0, x0[1]], f([x0[0] + 1.0, x0[1]])),
        ([x0[0], x0[1] + 1.0], f([x0[0], x0[1] + 1.0])),
    ];
    let mut iterations = 0;
    while iterations < 4000 {
        simplex.sort_by(|p, q| p.1.total_cmp(&q.1));
        let (b, w) = (simplex[0], simplex[2]);
        let span = simplex[1..]
            .iter()
            .map(|p| (p.0[0] - b.0[0]).abs().max((p.0[1] - b.0[1]).abs()))
            .fold(0.0, f64::max);
        let fspread = (w.1 - b.1).abs();
        if span < 1e-7 && fspread <= 1e-13 * b.1.abs().max(1e-300) {
            break;
        }
        iterations += 1;
        let c = [
            0.5 * (simplex[0].0[0] + simplex[1].0[0]),
            0.5 * (simplex[0].0[1] + simplex[1].0[1]),
        ];
        let at = |t: f64| [c[0] + t * (w.0[0] - c[0]), c[1] + t * (w.0[1] - c[1])];
        let xr = at(-1.0);
        let fr = f(xr);
        if fr < simplex[0].1 {
            let xe = at(-2.0);
            let fe = f(xe);
            simplex[2] = if fe < fr { (xe, fe) } else { (xr, fr) };
        } else if fr < simplex[1].1 {
            simplex[2] = (xr, fr);
        } else {
            let (xc, fc) = if fr < w.1 {
                let x = at(-0.5);
                (x, f(x))
            } else {
                let x = at(0.5);
                (x, f(x))
            };
            if fc < fr.min(w.1) {
                simplex[2] = (xc, fc);
            } else {
                for p in simplex.iter_mut().skip(1) {
                    let x = [
                        b.0[0] + 0.5 * (p.0[0] - b.0[0]),
                        b.0[1] + 0.5 * (p.0[1] - b.0[1]),
                    ];
                    *p = (x, f(x));
                }
            }
        }
    }
    simplex.sort_by(|p, q| p.1.total_cmp(&q.1));
    let x = simplex[0].0;
    let (c3, vinf, theta) = transfer_eval(kernel, dep_body, arr_body, x[0], x[1], ttype, mu_sun)
        .ok_or("the optimum left the feasible region")?;
    Ok(TransferOptimum {
        dep_jd_tdb: x[0],
        arr_jd_tdb: x[1],
        c3_km2s2: c3,
        vinf_arr_kms: vinf,
        transfer_angle_deg: theta,
        iterations,
    })
}

/// Closed-form departure C3 (m²/s²) of the Hohmann transfer between two coplanar circular
/// orbits `r1 → r2` about `mu`: `C3 = (v_transfer,peri − v_circ,1)²`. This is the global
/// minimum-energy departure for circular-to-circular transfer — the porkchop floor.
pub fn hohmann_departure_c3(r1: f64, r2: f64, mu: f64) -> f64 {
    let a_t = 0.5 * (r1 + r2);
    let v_circ1 = (mu / r1).sqrt();
    let v_peri = (mu * (2.0 / r1 - 1.0 / a_t)).sqrt();
    let vinf = v_peri - v_circ1;
    vinf * vinf
}

/// Time of flight (s) of the Hohmann transfer `r1 → r2` about `mu`: half the transfer-ellipse
/// period, `π·√(a_t³/μ)` with `a_t = (r1+r2)/2`.
pub fn hohmann_tof(r1: f64, r2: f64, mu: f64) -> f64 {
    let a_t = 0.5 * (r1 + r2);
    PI * (a_t * a_t * a_t / mu).sqrt()
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::f64::consts::PI;

    fn approx(a: f64, b: f64, tol: f64) {
        assert!((a - b).abs() <= tol, "expected {a} ≈ {b} (tol {tol})");
    }

    // ----- 1. Impulsive maneuver + covariance -----

    #[test]
    fn impulsive_eci_dv_jumps_velocity_and_adds_only_velocity_covariance() {
        let state = [7.0e6, 0.0, 0.0, 0.0, 7.5e3, 0.0];
        let mut cov = [[0.0f64; 6]; 6];
        for (i, row) in cov.iter_mut().enumerate() {
            row[i] = (i as f64 + 1.0) * 10.0; // distinct diagonal
        }
        let man = ImpulsiveManeuver {
            dv: [10.0, -5.0, 2.0],
            frame: ManeuverFrame::Eci,
            exec_cov: [[0.04, 0.0, 0.0], [0.0, 0.09, 0.0], [0.0, 0.0, 0.01]],
        };
        let (out, p) = apply_impulse(state, cov, &man);
        // Position unchanged, velocity jumps by ΔV exactly.
        approx(out[0], 7.0e6, 0.0);
        approx(out[3], 10.0, 1e-12);
        approx(out[4], 7.5e3 - 5.0, 1e-9);
        approx(out[5], 2.0, 1e-12);
        // Position block of P unchanged; velocity block gains exactly exec_cov.
        for i in 0..3 {
            for j in 0..3 {
                approx(p[i][j], cov[i][j], 1e-12);
            }
        }
        approx(p[3][3], cov[3][3] + 0.04, 1e-12);
        approx(p[4][4], cov[4][4] + 0.09, 1e-12);
        approx(p[5][5], cov[5][5] + 0.01, 1e-12);
    }

    #[test]
    fn lvlh_execution_covariance_is_rotated_trace_preserving_and_symmetric() {
        // A non-trivial inclined state so the LVLH frame is genuinely rotated.
        let state = [7.0e6, 1.0e6, 5.0e5, 1.0e3, 7.2e3, 1.5e3];
        let cov = [[0.0f64; 6]; 6];
        let (a, b, c) = (0.05, 0.20, 0.02);
        let man = ImpulsiveManeuver {
            dv: [3.0, 0.0, 0.0],
            frame: ManeuverFrame::Lvlh,
            exec_cov: [[a, 0.0, 0.0], [0.0, b, 0.0], [0.0, 0.0, c]],
        };
        let (_out, p) = apply_impulse(state, cov, &man);
        // Rotation preserves the trace of the covariance added to the velocity block.
        let tr = p[3][3] + p[4][4] + p[5][5];
        approx(tr, a + b + c, 1e-12);
        // The added velocity-block covariance is symmetric.
        approx(p[3][4], p[4][3], 1e-15);
        approx(p[3][5], p[5][3], 1e-15);
        approx(p[4][5], p[5][4], 1e-15);
        // And the LVLH→ECI columns are orthonormal.
        let r = [state[0], state[1], state[2]];
        let v = [state[3], state[4], state[5]];
        let m = lvlh_to_eci(r, v);
        for col in &m {
            approx(norm(*col), 1.0, 1e-12);
        }
        approx(dot(m[0], m[1]), 0.0, 1e-12);
        approx(dot(m[0], m[2]), 0.0, 1e-12);
        approx(dot(m[1], m[2]), 0.0, 1e-12);
    }

    // ----- 2. Finite burn vs Tsiolkovsky -----

    #[test]
    fn finite_burn_matches_tsiolkovsky_to_better_than_a_hundredth_percent() {
        // Hand check: ṁ = 500/(300·9.80665) = 0.169970 kg/s; mf = 1000 − 16.9970 = 983.003 kg;
        // Δv = 300·9.80665·ln(1000/983.003) = 2941.995·0.0171419 ≈ 50.43 m/s.
        let burn = FiniteBurn {
            thrust_n: 500.0,
            isp_s: 300.0,
            m0_kg: 1000.0,
            burn_s: 100.0,
            dir: [0.0, 1.0, 0.0],
        };
        let res = integrate_finite_burn(&burn, 2000).unwrap();
        let mdot = 500.0 / (300.0 * G0);
        approx(res.mf_kg, 1000.0 - mdot * 100.0, 1e-9);
        approx(
            res.tsiolkovsky_ms,
            tsiolkovsky(300.0, 1000.0, res.mf_kg),
            0.0,
        );
        // Closed-form sanity on the magnitude itself.
        approx(res.tsiolkovsky_ms, 50.43, 0.05);
        // RK4 of v̇ = T/m reproduces the exact rocket equation to well under 0.01 %.
        assert!(
            res.rel_err < 1e-4,
            "finite-burn ΔV rel err {} ≥ 1e-4",
            res.rel_err
        );
    }

    #[test]
    fn finite_burn_errors_when_mass_is_exhausted() {
        let burn = FiniteBurn {
            thrust_n: 5000.0,
            isp_s: 200.0,
            m0_kg: 100.0,
            burn_s: 100.0,
            dir: [1.0, 0.0, 0.0],
        };
        assert!(integrate_finite_burn(&burn, 100).is_err());
    }

    // ----- 3. Kepler propagator self-check -----

    #[test]
    fn kepler_universal_round_trips_forward_then_back() {
        let r0 = [7.0e6, 0.0, 0.0];
        let v0 = [0.0, 8.0e3, 1.0e3];
        let mu = crate::orbit::MU_EARTH;
        let (r1, v1) = kepler_universal(r0, v0, 2000.0, mu);
        let (r2, v2) = kepler_universal(r1, v1, -2000.0, mu);
        approx(norm(sub(r2, r0)), 0.0, 1e-3);
        approx(norm(sub(v2, v0)), 0.0, 1e-6);
    }

    // ----- 4. Lambert: round-trip against exact two-body truth -----

    #[test]
    fn lambert_recovers_the_velocities_of_a_known_two_body_arc() {
        let mu = crate::orbit::MU_EARTH;
        let r1 = [7.0e6, 0.0, 0.0];
        let v1_true = [0.0, 8.0e3, 1.0e3]; // h = r×v has +z ⇒ prograde
        let tof = 2000.0;
        // Exact two-body truth: where does this arc end, and at what velocity?
        let (r2, v2_true) = kepler_universal(r1, v1_true, tof, mu);
        // Lambert must reconstruct both boundary velocities.
        let (v1, v2) = lambert(r1, r2, tof, mu, true).unwrap();
        approx(norm(sub(v1, v1_true)), 0.0, 1e-3);
        approx(norm(sub(v2, v2_true)), 0.0, 1e-3);
    }

    #[test]
    fn lambert_output_propagates_back_onto_r2() {
        // Independent of the orbit fixture: any Lambert solution, when propagated, must hit r2.
        let mu = MU_SUN;
        let r1 = [1.0 * AU_M, 0.2 * AU_M, 0.0];
        let r2 = [-0.3 * AU_M, 1.3 * AU_M, 0.05 * AU_M];
        let tof = 200.0 * 86400.0;
        let (v1, v2) = lambert(r1, r2, tof, mu, true).unwrap();
        let (r_end, v_end) = kepler_universal(r1, v1, tof, mu);
        approx(norm(sub(r_end, r2)), 0.0, 1.0e3); // < 1 km on a ~1.5e11 m arc
        approx(norm(sub(v_end, v2)), 0.0, 1e-3);
    }

    #[test]
    fn lambert_rejects_collinear_endpoints() {
        let mu = crate::orbit::MU_EARTH;
        let r1 = [7.0e6, 0.0, 0.0];
        let r2 = [-9.0e6, 0.0, 0.0]; // exactly 180°, plane undefined
        assert!(lambert(r1, r2, 3000.0, mu, true).is_err());
    }

    // ----- 5. Porkchop: round-trip consistency + Hohmann floor -----

    #[test]
    fn porkchop_cells_are_round_trip_consistent_and_min_is_near_the_hohmann_floor() {
        // Synthetic coplanar-circular heliocentric system: 1.0 AU → 1.524 AU.
        let r1 = AU_M;
        let r2 = 1.524 * AU_M;
        let tof_h = hohmann_tof(r1, r2, MU_SUN);
        let dep = CircularCoplanarBody {
            radius_m: r1,
            phase0_rad: 0.0,
            mu_central: MU_SUN,
        };
        // Place the target so a 0 → tof_h transfer sweeps ~180° (the Hohmann geometry),
        // but we never sample exactly there (the guard / grid offsets keep us off the singularity).
        let n2 = (MU_SUN / r2.ppowi(3)).sqrt();
        let arr = CircularCoplanarBody {
            radius_m: r2,
            phase0_rad: PI - n2 * tof_h,
            mu_central: MU_SUN,
        };
        let day = 86400.0;
        let dep_epochs: Vec<f64> = (-20..=20).map(|k| k as f64 * 5.0 * day).collect();
        let arr_epochs: Vec<f64> = (0..=40)
            .map(|k| tof_h - 100.0 * day + k as f64 * 5.0 * day)
            .collect();
        let grid = porkchop(&dep, &arr, &dep_epochs, &arr_epochs, MU_SUN);

        // Every finite cell's departure velocity must propagate two-body back onto r2.
        let mut finite_cells = 0;
        for (i, &td) in dep_epochs.iter().enumerate() {
            for (j, &ta) in arr_epochs.iter().enumerate() {
                if grid.c3_km2s2[i][j].is_finite() {
                    finite_cells += 1;
                    let tof = ta - td;
                    let rr1 = dep.position(td);
                    let rr2 = arr.position(ta);
                    let (v1, _) = lambert(rr1, rr2, tof, MU_SUN, true).unwrap();
                    let (r_end, _) = kepler_universal(rr1, v1, tof, MU_SUN);
                    assert!(
                        norm(sub(r_end, rr2)) < 5.0e3,
                        "cell ({i},{j}) round-trip miss {} m",
                        norm(sub(r_end, rr2))
                    );
                }
            }
        }
        assert!(finite_cells > 100, "too few finite cells: {finite_cells}");

        // The closed-form Hohmann departure C3 is the theoretical floor.
        let c3_floor = hohmann_departure_c3(r1, r2, MU_SUN) / 1.0e6; // km²/s²
        let (min_c3, mi, mj) = grid.min_c3().unwrap();
        // Floor property: no cell beats Hohmann (small numerical slack).
        assert!(
            min_c3 >= c3_floor - 1e-6,
            "grid min C3 {min_c3} below Hohmann floor {c3_floor}"
        );
        // The finite grid samples near (but not at) 180°, so the minimum is just above the floor.
        assert!(
            min_c3 < c3_floor * 1.05,
            "grid min C3 {min_c3} not within 5% of Hohmann floor {c3_floor}"
        );
        // The minimizing cell's TOF is near the Hohmann TOF.
        let tof_min = arr_epochs[mj] - dep_epochs[mi];
        approx(tof_min, tof_h, 60.0 * day);
    }

    #[test]
    fn porkchop_json_round_trips_with_matching_dimensions() {
        let dep = CircularCoplanarBody {
            radius_m: AU_M,
            phase0_rad: 0.0,
            mu_central: MU_SUN,
        };
        let arr = CircularCoplanarBody {
            radius_m: 1.524 * AU_M,
            phase0_rad: 1.0,
            mu_central: MU_SUN,
        };
        let day = 86400.0;
        let dep_epochs: Vec<f64> = (0..5).map(|k| k as f64 * 10.0 * day).collect();
        let arr_epochs: Vec<f64> = (0..6).map(|k| (200 + k * 10) as f64 * day).collect();
        let grid = porkchop(&dep, &arr, &dep_epochs, &arr_epochs, MU_SUN);
        let json = grid.to_json();
        let parsed: serde_json::Value = serde_json::from_str(&json).unwrap();
        assert_eq!(parsed["dep_epochs_s"].as_array().unwrap().len(), 5);
        assert_eq!(parsed["arr_epochs_s"].as_array().unwrap().len(), 6);
        assert_eq!(parsed["c3_km2s2"].as_array().unwrap().len(), 5);
        assert_eq!(parsed["c3_km2s2"][0].as_array().unwrap().len(), 6);
        assert_eq!(parsed["vinf_arr_kms"].as_array().unwrap().len(), 5);
    }
}