kshana 0.36.0

// SPDX-License-Identifier: AGPL-3.0-only
//! Hyperbolic-flyby **B-plane targeting** and **patched-conic gravity assist**.
//!
//! A flyby is a hyperbolic encounter in the body-centred frame. Its geometry is fixed by
//! the gravitational parameter `μ`, the hyperbolic-excess speed `v∞` (the speed far from
//! the body, where the trajectory is asymptotic), and the periapsis radius `r_p`:
//!
//! ```text
//!   a = −μ / v∞²                 (semi-major axis, negative for a hyperbola)
//!   e = 1 + r_p · v∞² / μ        (eccentricity, > 1)
//!   δ = 2·asin(1/e)              (turn / deflection angle)
//!   |B| = |a|·√(e²−1) = r_p·√((e+1)/(e−1))   (impact parameter = semi-minor axis)
//! ```
//!
//! The **B-plane** is the plane through the body centre perpendicular to the incoming
//! asymptote `Ŝ`. With a reference pole `p̂`, the in-plane axes are `T̂ = Ŝ × p̂` (unit)
//! and `R̂ = Ŝ × T̂`; the aim point is given by the two scalars `B·T̂`, `B·R̂`, and
//! `|B|² = (B·T̂)² + (B·R̂)²`.
//!
//! In a **gravity assist**, the body-relative speed `|v∞|` is conserved while its
//! direction rotates by `δ` in the flyby plane. In the heliocentric frame the spacecraft
//! velocity is `v = v_planet + v∞`, so the assist imparts
//! `Δv = v∞_out − v∞_in` with magnitude `2·v∞·sin(δ/2)` (≤ `2·v∞`, the head-on limit) at
//! **no propellant cost**. The energy is borrowed from the planet's heliocentric motion;
//! the conserved bookkeeping quantity is the **Tisserand parameter** with respect to the
//! planet's (circular) orbit of radius `a_P`,
//!
//! ```text
//!   T_P = a_P/a + 2·√((a/a_P)(1−e²))·cos i,        with   v∞² = v_circ²·(3 − T_P),
//! ```
//!
//! which is invariant across the encounter because `|v∞|` is.
//!
//! Scope (honest): patched-conic two-body flyby on a circular planetary orbit — no
//! finite-sphere-of-influence transition modelling, no third-body perturbations during
//! the encounter, and no ephemeris. The unit tests check the closed-form flyby identities,
//! the B-plane decomposition and Tisserand invariance across a v∞-preserving deflection;
//! `tests/bplane_gmat_oracle.rs` compares [`flyby_from_state`] (B·T̂, B·R̂, |B|, the turn
//! angle and both asymptotes) with the General Mission Analysis Tool (GMAT) R2026a on 192
//! hyperbolic states, and `tests/bplane_heliocentric_oracle.rs` compares [`assist_delta_v`],
//! [`elements_aei`] and [`tisserand`] with GMAT, sbpy and published Tisserand values.
//!
//! References:
//! - R. H. Battin, *An Introduction to the Mathematics and Methods of Astrodynamics*,
//!   rev. ed., AIAA 1999, §6 (hyperbolic orbits, the B-plane).
//! - D. A. Vallado, *Fundamentals of Astrodynamics and Applications*, 4th ed., §12
//!   (interplanetary trajectories, gravity assist, the Tisserand criterion).

use crate::frames::Vec3;
use crate::portable_math::PortableFloat;

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]
}
fn sub(a: Vec3, b: Vec3) -> Vec3 {
    [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
}

/// Hyperbolic semi-major axis `a = −μ/v∞²` (m; negative). `v_inf` is the
/// hyperbolic-excess speed (m/s).
///
/// Requires `v_inf != 0`: a zero-energy encounter is not a hyperbola and the
/// closed form returns `−∞` rather than an error. The signature is kept
/// infallible because the whole module is a set of closed forms over an
/// already-hyperbolic encounter; the precondition is checked in debug builds.
pub fn hyperbolic_sma(mu: f64, v_inf: f64) -> f64 {
    debug_assert!(v_inf != 0.0, "hyperbolic_sma needs a non-zero v_inf");
    -mu / (v_inf * v_inf)
}

/// Flyby eccentricity `e = 1 + r_p·v∞²/μ` (> 1) for periapsis radius `r_p` (m).
pub fn flyby_eccentricity(mu: f64, v_inf: f64, r_p: f64) -> f64 {
    1.0 + r_p * v_inf * v_inf / mu
}

/// Turn (deflection) angle `δ = 2·asin(1/e)` (rad) for flyby eccentricity `e`.
///
/// Requires `e >= 1` (the flyby must be hyperbolic, or parabolic at the limit,
/// where `turn_angle(1.0) = π`). For `e < 1` the argument of `asin` leaves
/// `[-1, 1]` and the result is `NaN`, which then survives every arithmetic
/// operation downstream — so the precondition is checked in debug builds.
pub fn turn_angle(e: f64) -> f64 {
    debug_assert!(e >= 1.0, "turn_angle needs a hyperbolic eccentricity");
    2.0 * (1.0 / e).pasin()
}

/// Impact parameter (B-plane aim radius) `|B| = |a|·√(e²−1)` (m).
///
/// Requires `v_inf != 0` (see [`hyperbolic_sma`]): at `v_inf = 0` the two
/// factors are `∞` and `0` and the product is `NaN`.
pub fn impact_parameter(mu: f64, v_inf: f64, r_p: f64) -> f64 {
    let a = hyperbolic_sma(mu, v_inf);
    let e = flyby_eccentricity(mu, v_inf, r_p);
    a.abs() * (e * e - 1.0).sqrt()
}

/// Rotate `v` by angle `delta` about the unit axis `axis` (Rodrigues' rotation). Used to
/// deflect the incoming `v∞` into the outgoing `v∞` across a flyby.
pub fn deflect(v: Vec3, delta: f64, axis: Vec3) -> Vec3 {
    let k = unit(axis);
    let (s, c) = (delta.psin(), delta.pcos());
    let kv = dot(k, v);
    // v·cosδ + (k×v)·sinδ + k·(k·v)·(1−cosδ)
    let kxv = cross(k, v);
    [
        v[0] * c + kxv[0] * s + k[0] * kv * (1.0 - c),
        v[1] * c + kxv[1] * s + k[1] * kv * (1.0 - c),
        v[2] * c + kxv[2] * s + k[2] * kv * (1.0 - c),
    ]
}

/// The heliocentric velocity change imparted by a flyby: `Δv = v∞_out − v∞_in` (the
/// planet velocity cancels). Its magnitude is `2·v∞·sin(δ/2)`.
pub fn assist_delta_v(v_inf_in: Vec3, v_inf_out: Vec3) -> Vec3 {
    sub(v_inf_out, v_inf_in)
}

/// Construct the B-plane in-plane axes `(T̂, R̂)` from the incoming-asymptote unit
/// vector `s_hat` and a reference pole `pole`: `T̂ = Ŝ × p̂` (unit), `R̂ = Ŝ × T̂`.
pub fn bplane_frame(s_hat: Vec3, pole: Vec3) -> (Vec3, Vec3) {
    let t = unit(cross(s_hat, pole));
    let r = cross(s_hat, t);
    (t, r)
}

/// Decompose a B-vector onto the B-plane axes: `(B·T̂, B·R̂)`.
pub fn bplane_components(b_vec: Vec3, t_hat: Vec3, r_hat: Vec3) -> (f64, f64) {
    (dot(b_vec, t_hat), dot(b_vec, r_hat))
}

/// Heliocentric osculating elements `(a, e, i)` from a state `(r, v)` and `μ`:
/// `a` (vis-viva, m), `e` (eccentricity-vector magnitude), `i` (rad, from the reference
/// plane). Two-body; valid for a bound (`a > 0`) heliocentric arc.
pub fn elements_aei(r: Vec3, v: Vec3, mu: f64) -> (f64, f64, f64) {
    let rm = norm(r);
    let v2 = dot(v, v);
    let a = 1.0 / (2.0 / rm - v2 / mu);
    let h = cross(r, v);
    let hm = norm(h);
    let i = (h[2] / hm).pacos();
    // eccentricity vector e = ((v²−μ/r)·r − (r·v)·v)/μ
    let rv = dot(r, v);
    let e_vec = [
        ((v2 - mu / rm) * r[0] - rv * v[0]) / mu,
        ((v2 - mu / rm) * r[1] - rv * v[1]) / mu,
        ((v2 - mu / rm) * r[2] - rv * v[2]) / mu,
    ];
    (a, norm(e_vec), i)
}

/// Tisserand parameter with respect to a planet on a circular orbit of radius `a_planet`:
/// `T_P = a_P/a + 2·√((a/a_P)(1−e²))·cos i`. Conserved across a patched-conic flyby
/// because it is a function of the (invariant) `v∞` alone.
pub fn tisserand(a_sc: f64, e_sc: f64, i_sc: f64, a_planet: f64) -> f64 {
    a_planet / a_sc + 2.0 * (a_sc / a_planet * (1.0 - e_sc * e_sc)).sqrt() * i_sc.pcos()
}

/// The B-plane description of a hyperbolic body-centred state, from [`flyby_from_state`].
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct FlybyState {
    /// Hyperbolic-excess speed `v∞ = √(v² − 2μ/r)` (state units).
    pub v_inf: f64,
    /// Periapsis radius `r_p = |a|·(e − 1)`.
    pub r_p: f64,
    /// Eccentricity (> 1).
    pub e: f64,
    /// Turn angle `δ = 2·asin(1/e)` (rad).
    pub turn_angle: f64,
    /// Impact parameter `|B|`.
    pub b_mag: f64,
    /// `B·T̂`.
    pub b_dot_t: f64,
    /// `B·R̂`.
    pub b_dot_r: f64,
    /// Incoming-asymptote unit vector `Ŝ` (the arrival direction of `v∞`).
    pub s_in: Vec3,
    /// Outgoing-asymptote unit vector: `Ŝ` deflected by `δ` about the orbit normal.
    pub s_out: Vec3,
    /// Unit orbit normal `ĥ = r × v / |r × v|`.
    pub h_hat: Vec3,
}

/// B-plane coordinates of a hyperbolic body-centred state `(r, v)` with gravitational
/// parameter `mu`, with `T̂ = Ŝ × p̂` built from the reference pole `pole` (the `z` axis of
/// the frame for the usual convention).
///
/// The incoming asymptote is `Ŝ = ê/e + √(1 − 1/e²)·(ĥ × ê)` (the arrival direction of the
/// hyperbolic-excess velocity), the B-vector is `|B|·(Ŝ × ĥ)` with `|B|` from
/// [`impact_parameter`], and the outgoing asymptote is `Ŝ` rotated by the turn angle about
/// `ĥ` with [`deflect`]. Units follow the inputs (metres and m³/s², or kilometres and
/// km³/s²). Returns `None` for a non-hyperbolic state (`e ≤ 1` or non-positive energy) or a
/// degenerate one (rectilinear, or `Ŝ` parallel to the pole).
pub fn flyby_from_state(r: Vec3, v: Vec3, mu: f64, pole: Vec3) -> Option<FlybyState> {
    let rm = norm(r);
    let v2 = dot(v, v);
    let energy2 = v2 - 2.0 * mu / rm;
    if energy2.is_nan() || energy2 <= 0.0 {
        return None;
    }
    let v_inf = energy2.sqrt();
    let h = cross(r, v);
    let hm = norm(h);
    if hm.is_nan() || hm <= 0.0 {
        return None;
    }
    let h_hat = [h[0] / hm, h[1] / hm, h[2] / hm];
    let rv = dot(r, v);
    let e_vec = [
        ((v2 - mu / rm) * r[0] - rv * v[0]) / mu,
        ((v2 - mu / rm) * r[1] - rv * v[1]) / mu,
        ((v2 - mu / rm) * r[2] - rv * v[2]) / mu,
    ];
    let e = norm(e_vec);
    if e.is_nan() || e <= 1.0 {
        return None;
    }
    let e_hat = [e_vec[0] / e, e_vec[1] / e, e_vec[2] / e];
    let n_hat = cross(h_hat, e_hat);
    let a = hyperbolic_sma(mu, v_inf);
    let r_p = a.abs() * (e - 1.0);
    let c = 1.0 / e;
    let s = (1.0 - c * c).sqrt();
    let s_in = [
        c * e_hat[0] + s * n_hat[0],
        c * e_hat[1] + s * n_hat[1],
        c * e_hat[2] + s * n_hat[2],
    ];
    if norm(cross(s_in, pole)) <= 1e-12 {
        return None;
    }
    let b_mag = impact_parameter(mu, v_inf, r_p);
    let b_dir = unit(cross(s_in, h_hat));
    let b_vec = [b_mag * b_dir[0], b_mag * b_dir[1], b_mag * b_dir[2]];
    let (t_hat, r_hat) = bplane_frame(s_in, pole);
    let (b_dot_t, b_dot_r) = bplane_components(b_vec, t_hat, r_hat);
    let delta = turn_angle(flyby_eccentricity(mu, v_inf, r_p));
    let s_out = deflect(s_in, delta, h_hat);
    Some(FlybyState {
        v_inf,
        r_p,
        e,
        turn_angle: delta,
        b_mag,
        b_dot_t,
        b_dot_r,
        s_in,
        s_out,
        h_hat,
    })
}

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

    fn approx(a: f64, b: f64, tol: f64) -> bool {
        (a - b).abs() <= tol
    }

    // Earth-ish μ (m³/s²) for the flyby-scalar tests.
    const MU: f64 = 3.986_004_418e14;

    #[test]
    fn flyby_from_state_recovers_the_closed_forms_at_periapsis() {
        // Periapsis state of a known flyby in the xy-plane, prograde about +z.
        let (v_inf, r_p) = (3000.0, 7.0e6);
        let v_p = (v_inf * v_inf + 2.0 * MU / r_p).sqrt();
        let pole = [0.3, 0.0, 1.0];
        let f = flyby_from_state([r_p, 0.0, 0.0], [0.0, v_p, 0.0], MU, pole).expect("hyperbolic");
        assert!(approx(f.v_inf, v_inf, 1e-6) && approx(f.r_p, r_p, 1e-3));
        let e = flyby_eccentricity(MU, v_inf, r_p);
        assert!(approx(f.e, e, 1e-12));
        assert!(approx(f.turn_angle, turn_angle(e), 1e-12));
        assert!(approx(f.b_mag, impact_parameter(MU, v_inf, r_p), 1e-3));
        assert!(approx(f.b_dot_t.phypot(f.b_dot_r), f.b_mag, 1e-3));
        let ang = dot(f.s_in, f.s_out).clamp(-1.0, 1.0).pacos();
        assert!(approx(ang, f.turn_angle, 1e-9));
        // The outgoing asymptote is the incoming one reflected through the apse line.
        let e_hat = [1.0, 0.0, 0.0];
        let refl = 2.0 * dot(f.s_in, e_hat);
        let mirrored = [refl - f.s_in[0], -f.s_in[1], -f.s_in[2]];
        assert!(norm(sub(mirrored, [-f.s_out[0], -f.s_out[1], -f.s_out[2]])) < 1e-12);
        // An elliptic state is rejected.
        let ell = flyby_from_state([r_p, 0.0, 0.0], [0.0, 7000.0, 0.0], MU, pole);
        assert!(ell.is_none());
    }

    #[test]
    fn flyby_scalars_satisfy_the_closed_forms() {
        let v_inf = 3000.0; // m/s
        let r_p = 7.0e6; // m
        let a = hyperbolic_sma(MU, v_inf);
        assert!(a < 0.0, "hyperbolic a must be negative");
        assert!(approx(a, -MU / (v_inf * v_inf), 1e-3));
        let e = flyby_eccentricity(MU, v_inf, r_p);
        assert!(e > 1.0, "flyby must be hyperbolic, e={e}");
        // turn angle identity sin(δ/2) = 1/e.
        let delta = turn_angle(e);
        assert!(approx((delta / 2.0).psin(), 1.0 / e, 1e-12));
        // impact parameter: two closed forms agree.
        let b1 = impact_parameter(MU, v_inf, r_p);
        let b2 = r_p * ((e + 1.0) / (e - 1.0)).sqrt();
        assert!(approx(b1, b2, 1e-3), "|B|: {b1} vs {b2}");
        assert!(b1 > r_p, "gravitational focusing ⇒ |B| > r_p");
    }

    #[test]
    fn turn_angle_decreases_with_periapsis_radius() {
        let v_inf = 2500.0;
        let close = turn_angle(flyby_eccentricity(MU, v_inf, 6.6e6));
        let far = turn_angle(flyby_eccentricity(MU, v_inf, 5.0e7));
        assert!(
            close > far,
            "closer flyby must deflect more: {close} vs {far}"
        );
        // limits: r_p → ∞ ⇒ δ → 0; e → 1⁺ ⇒ δ → π.
        assert!(turn_angle(flyby_eccentricity(MU, v_inf, 1.0e12)) < 1e-3);
        assert!(turn_angle(1.0 + 1e-9) > PI - 1e-3);
    }

    #[test]
    fn deflection_preserves_speed_and_rotates_by_delta() {
        let v_in: Vec3 = [3000.0, 500.0, -200.0];
        let axis: Vec3 = [0.2, -0.4, 1.0];
        // axis component along v is irrelevant; use the perpendicular plane normal.
        let normal = unit(cross(v_in, axis));
        let delta = 0.7;
        let v_out = deflect(v_in, delta, normal);
        assert!(approx(norm(v_out), norm(v_in), 1e-6), "speed not preserved");
        let ang = (dot(v_in, v_out) / (norm(v_in) * norm(v_out))).pacos();
        assert!(approx(ang, delta, 1e-9), "rotation angle {ang} vs {delta}");
    }

    #[test]
    fn assist_delta_v_magnitude_is_two_vinf_sin_half_delta() {
        let v_inf = 4000.0;
        let v_in: Vec3 = [v_inf, 0.0, 0.0];
        let normal: Vec3 = [0.0, 0.0, 1.0];
        for &delta in &[0.3_f64, 1.0, 2.0, PI] {
            let v_out = deflect(v_in, delta, normal);
            let dv = norm(assist_delta_v(v_in, v_out));
            assert!(
                approx(dv, 2.0 * v_inf * (delta / 2.0).psin(), 1e-6),
                "Δv {dv} vs {}",
                2.0 * v_inf * (delta / 2.0).psin()
            );
            assert!(dv <= 2.0 * v_inf + 1e-6, "Δv exceeds the head-on limit");
        }
    }

    #[test]
    fn bplane_frame_is_orthonormal_and_decomposes_b() {
        let s_hat = unit([1.0, 0.3, -0.5]);
        let pole: Vec3 = [0.0, 0.0, 1.0];
        let (t, r) = bplane_frame(s_hat, pole);
        // orthonormal and both perpendicular to S.
        assert!(approx(norm(t), 1.0, 1e-12) && approx(norm(r), 1.0, 1e-12));
        assert!(approx(dot(t, r), 0.0, 1e-12));
        assert!(approx(dot(t, s_hat), 0.0, 1e-12) && approx(dot(r, s_hat), 0.0, 1e-12));
        // a B-vector in the plane decomposes with |B|² = (B·T)² + (B·R)².
        let b_vec = [
            3.0 * t[0] + 4.0 * r[0],
            3.0 * t[1] + 4.0 * r[1],
            3.0 * t[2] + 4.0 * r[2],
        ];
        let (bt, br) = bplane_components(b_vec, t, r);
        assert!(approx(bt, 3.0, 1e-9) && approx(br, 4.0, 1e-9));
        assert!(approx(norm(b_vec) * norm(b_vec), bt * bt + br * br, 1e-9));
    }

    #[test]
    fn tisserand_is_invariant_across_a_vinf_preserving_deflection() {
        // Normalised system: μ_sun = 1, planet circular orbit a_P = 1 ⇒ v_circ = 1.
        let mu_sun = 1.0_f64;
        let a_p = 1.0_f64;
        let v_circ = (mu_sun / a_p).sqrt();
        let r: Vec3 = [a_p, 0.0, 0.0];
        let v_planet: Vec3 = [0.0, v_circ, 0.0];
        let v_inf_in: Vec3 = [0.10, 0.20, 0.05];
        let v_inf_mag = norm(v_inf_in);

        let v_helio_in = [
            v_planet[0] + v_inf_in[0],
            v_planet[1] + v_inf_in[1],
            v_planet[2] + v_inf_in[2],
        ];
        let (a0, e0, i0) = elements_aei(r, v_helio_in, mu_sun);
        assert!(a0 > 0.0, "spacecraft must be on a bound heliocentric orbit");
        let t0 = tisserand(a0, e0, i0, a_p);

        // Deflect v∞ by an arbitrary angle about an arbitrary axis (speed preserved).
        let axis = unit([0.3, 1.0, -0.2]);
        let v_inf_out = deflect(v_inf_in, 1.1, axis);
        let v_helio_out = [
            v_planet[0] + v_inf_out[0],
            v_planet[1] + v_inf_out[1],
            v_planet[2] + v_inf_out[2],
        ];
        let (a1, e1, i1) = elements_aei(r, v_helio_out, mu_sun);
        let t1 = tisserand(a1, e1, i1, a_p);

        // Invariance: the assist changes a, e, i but not the Tisserand parameter.
        assert!(
            approx(t0, t1, 1e-9),
            "Tisserand not conserved: {t0} vs {t1}"
        );
        // And the v∞ link: T_P = 3 − (v∞/v_circ)².
        let t_from_vinf = 3.0 - (v_inf_mag / v_circ).ppowi(2);
        assert!(
            approx(t0, t_from_vinf, 1e-9),
            "T_P {t0} vs 3−(v∞/v_c)² {t_from_vinf}"
        );
        // sanity: the orbit really did change.
        assert!((a0 - a1).abs() + (e0 - e1).abs() + (i0 - i1).abs() > 1e-3);
    }
}