ballistics_engine/

trajectory_integration.rs

//! Advanced trajectory integration methods (RK4, RK45)
//!
//! This module provides production-grade numerical integration for ballistic trajectories:
//! - RK4: 4th-order Runge-Kutta (fixed step)
//! - RK45: Dormand-Prince adaptive method (same as scipy.integrate.solve_ivp)
//!
//! MBA-155: Upstreamed from ballistics_rust for shared use

use crate::gravity::Gravity;
use nalgebra::{Vector3, Vector6};
use std::collections::HashMap;

use crate::derivatives::compute_derivatives_with_gravity;
use crate::wind::{WindSegment, WindSegmentError};
use crate::BallisticInputs;
use crate::DragModel;

const RK45_MIN_STEP: f64 = 1e-6;
const RK45_DEFAULT_TOLERANCE: f64 = 1e-6;
const RK45_SAFETY_FACTOR: f64 = 0.9;
const RK45_MIN_SCALE: f64 = 0.1;
const RK45_MAX_SCALE: f64 = 2.0;

#[derive(Clone, Copy)]
struct Rk45Control {
    tolerance: f64,
    min_step: f64,
    max_step: f64,
    max_trials: usize,
}

struct Rk45AcceptedStep {
    state: Vector6<f64>,
    used_dt: f64,
    next_dt: f64,
    error: f64,
    trials: usize,
}

fn wind_vector_for_range(range_m: f64, wind_segments: &[WindSegment]) -> Vector3<f64> {
    if range_m.is_nan() {
        return Vector3::zeros();
    }
    for seg in wind_segments {
        if range_m < seg.until_m {
            let wind_speed_mps = seg.speed_kmh * 0.2777778; // km/h to m/s
            let wind_angle_rad = seg.angle_deg.to_radians();
            // MBA-728: per-segment vertical passes straight through (not derived from
            // speed/angle), matching wind::WindSock::calc_vec.
            return crate::wind::wind_vector(wind_speed_mps, wind_angle_rad, seg.vertical_mps);
        }
    }
    Vector3::zeros()
}

/// RK4 integration step
fn rk4_step_with_gravity(
    state: &Vector6<f64>,
    t: f64,
    dt: f64,
    params: &TrajectoryParams,
    inputs: &BallisticInputs,
    gravity: Gravity,
) -> Vector6<f64> {
    // RK4 integration
    let k1 = compute_derivatives_vec_with_gravity(state, t, params, inputs, gravity);
    let k2 = compute_derivatives_vec_with_gravity(
        &(state + dt * 0.5 * k1),
        t + dt * 0.5,
        params,
        inputs,
        gravity,
    );
    let k3 = compute_derivatives_vec_with_gravity(
        &(state + dt * 0.5 * k2),
        t + dt * 0.5,
        params,
        inputs,
        gravity,
    );
    let k4 =
        compute_derivatives_vec_with_gravity(&(state + dt * k3), t + dt, params, inputs, gravity);

    state + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
}

/// Weighted RMS error for a mixed position/velocity state.
///
/// Each component is scaled independently so a large downrange position cannot hide an error in
/// a near-zero lateral velocity (and vice versa). The caller's tolerance therefore acts as both
/// an absolute and relative tolerance in each component's own unit.
pub(crate) fn rk45_error_norm(
    state: &Vector6<f64>,
    fifth_order: &Vector6<f64>,
    fourth_order: &Vector6<f64>,
) -> f64 {
    let scaled_error_squared: f64 = (0..6)
        .map(|index| {
            let scale = 1.0 + state[index].abs().max(fifth_order[index].abs());
            ((fifth_order[index] - fourth_order[index]) / scale).powi(2)
        })
        .sum();

    (scaled_error_squared / 6.0).sqrt()
}

/// Adaptive RK45 integration step (Dormand-Prince method)
#[allow(clippy::too_many_arguments)]
fn rk45_step_with_gravity(
    state: &Vector6<f64>,
    t: f64,
    dt: f64,
    params: &TrajectoryParams,
    inputs: &BallisticInputs,
    tol: f64,
    gravity: Gravity,
) -> (Vector6<f64>, f64, f64) {
    // Dormand-Prince coefficients (same as scipy.integrate.solve_ivp RK45)
    const A21: f64 = 1.0 / 5.0;
    const A31: f64 = 3.0 / 40.0;
    const A32: f64 = 9.0 / 40.0;
    const A41: f64 = 44.0 / 45.0;
    const A42: f64 = -56.0 / 15.0;
    const A43: f64 = 32.0 / 9.0;
    const A51: f64 = 19372.0 / 6561.0;
    const A52: f64 = -25360.0 / 2187.0;
    const A53: f64 = 64448.0 / 6561.0;
    const A54: f64 = -212.0 / 729.0;
    const A61: f64 = 9017.0 / 3168.0;
    const A62: f64 = -355.0 / 33.0;
    const A63: f64 = 46732.0 / 5247.0;
    const A64: f64 = 49.0 / 176.0;
    const A65: f64 = -5103.0 / 18656.0;
    const A71: f64 = 35.0 / 384.0;
    const A73: f64 = 500.0 / 1113.0;
    const A74: f64 = 125.0 / 192.0;
    const A75: f64 = -2187.0 / 6784.0;
    const A76: f64 = 11.0 / 84.0;

    // 5th order coefficients
    const B1: f64 = 35.0 / 384.0;
    const B3: f64 = 500.0 / 1113.0;
    const B4: f64 = 125.0 / 192.0;
    const B5: f64 = -2187.0 / 6784.0;
    const B6: f64 = 11.0 / 84.0;

    // 4th order coefficients (for error estimation)
    const B1_ERR: f64 = 5179.0 / 57600.0;
    const B3_ERR: f64 = 7571.0 / 16695.0;
    const B4_ERR: f64 = 393.0 / 640.0;
    const B5_ERR: f64 = -92097.0 / 339200.0;
    const B6_ERR: f64 = 187.0 / 2100.0;
    const B7_ERR: f64 = 1.0 / 40.0;

    // Compute stages
    let k1 = compute_derivatives_vec_with_gravity(state, t, params, inputs, gravity);
    let k2 = compute_derivatives_vec_with_gravity(
        &(state + dt * A21 * k1),
        t + dt * 0.2,
        params,
        inputs,
        gravity,
    );
    let k3 = compute_derivatives_vec_with_gravity(
        &(state + dt * (A31 * k1 + A32 * k2)),
        t + dt * 0.3,
        params,
        inputs,
        gravity,
    );
    let k4 = compute_derivatives_vec_with_gravity(
        &(state + dt * (A41 * k1 + A42 * k2 + A43 * k3)),
        t + dt * 0.8,
        params,
        inputs,
        gravity,
    );
    let k5 = compute_derivatives_vec_with_gravity(
        &(state + dt * (A51 * k1 + A52 * k2 + A53 * k3 + A54 * k4)),
        t + dt * 8.0 / 9.0,
        params,
        inputs,
        gravity,
    );
    let k6 = compute_derivatives_vec_with_gravity(
        &(state + dt * (A61 * k1 + A62 * k2 + A63 * k3 + A64 * k4 + A65 * k5)),
        t + dt,
        params,
        inputs,
        gravity,
    );
    let k7 = compute_derivatives_vec_with_gravity(
        &(state + dt * (A71 * k1 + A73 * k3 + A74 * k4 + A75 * k5 + A76 * k6)),
        t + dt,
        params,
        inputs,
        gravity,
    );

    // 5th order solution
    let y_new = state + dt * (B1 * k1 + B3 * k3 + B4 * k4 + B5 * k5 + B6 * k6);

    // 4th order solution for error estimate
    let y_err = state
        + dt * (B1_ERR * k1 + B3_ERR * k3 + B4_ERR * k4 + B5_ERR * k5 + B6_ERR * k6 + B7_ERR * k7);

    let error = rk45_error_norm(state, &y_new, &y_err);

    // Dormand-Prince 5(4) controls both accepted and rejected trials with a fifth-root scale.
    let step_scale = if !error.is_finite() || !tol.is_finite() || tol <= 0.0 {
        RK45_MIN_SCALE
    } else if error == 0.0 {
        RK45_MAX_SCALE
    } else {
        (RK45_SAFETY_FACTOR * (tol / error).powf(0.2)).clamp(RK45_MIN_SCALE, RK45_MAX_SCALE)
    };
    let dt_new = dt * step_scale;

    (y_new, dt_new, error)
}

/// Retry an RK45 step from the same state until its embedded error estimate is acceptable.
///
/// `Err(n)` means no finite acceptable candidate was found in `n` trials. Rejected candidates
/// never escape this function, so callers cannot accidentally advance state or time with one.
#[allow(clippy::too_many_arguments)]
fn adaptive_rk45_step_with_gravity(
    state: &Vector6<f64>,
    t: f64,
    initial_dt: f64,
    params: &TrajectoryParams,
    inputs: &BallisticInputs,
    control: Rk45Control,
    gravity: Gravity,
) -> Result<Rk45AcceptedStep, usize> {
    let mut trial_dt = initial_dt;

    for trials in 1..=control.max_trials {
        let (new_state, suggested_dt, error) = rk45_step_with_gravity(
            state,
            t,
            trial_dt,
            params,
            inputs,
            control.tolerance,
            gravity,
        );
        let candidate_is_finite = error.is_finite()
            && suggested_dt.is_finite()
            && new_state.iter().all(|value| value.is_finite());
        let next_dt = suggested_dt.min(control.max_step).max(control.min_step);

        if candidate_is_finite && (error <= control.tolerance || trial_dt <= control.min_step) {
            return Ok(Rk45AcceptedStep {
                state: new_state,
                used_dt: trial_dt,
                next_dt,
                error,
                trials,
            });
        }

        if trial_dt <= control.min_step {
            return Err(trials);
        }
        trial_dt = next_dt;
    }

    Err(control.max_trials)
}

/// Parameters for trajectory computation
pub struct TrajectoryParams {
    pub mass_kg: f64,
    pub bc: f64,
    pub drag_model: DragModel,
    /// Downrange wind zones, normalized by `until_distance_m` when integration begins.
    pub wind_segments: Vec<WindSegment>,
    /// Dual-mode atmosphere tuple consumed by `compute_derivatives`:
    /// **Standard** `(base_alt_m, base_temp_c, base_pressure_hPa, base_density_ratio)` — note
    /// slot 3 is a density RATIO, NOT humidity, even though it rides in the `humidity` field;
    /// or **Direct** `(air_density, speed_of_sound, 0.0, 0.0)` — slots 2 and 3 are zero
    /// sentinels. A pressure of 0 that is not the direct-mode sentinel disables drag.
    pub atmos_params: (f64, f64, f64, f64),
    /// Earth rotation in level downrange/up/lateral axes. The derivative kernel projects it into
    /// the inclined shot frame using `shooting_angle` before applying Coriolis acceleration.
    pub omega_vector: Option<Vector3<f64>>,
    pub enable_spin_drift: bool,
    pub enable_magnus: bool,
    pub enable_coriolis: bool,
    pub target_distance_m: f64, // Target horizontal distance in meters
    /// Boundary-layer wind shear. Its height above the ground is the shot-frame y (the ground
    /// runs parallel to the line of sight, as for `ground_threshold`) with y = 0 at the launch
    /// point, taken as 1.5 m above the ground. `fast_integrate_with_segments` instead measures y
    /// from its inputs' `muzzle_height` datum.
    pub enable_wind_shear: bool,
    pub wind_shear_model: String,
    /// Site elevation above sea level, kept for API compatibility: nothing in this kernel reads
    /// it, and boundary-layer shear is relative to the local ground.
    pub shooter_altitude_m: f64,
    pub is_twist_right: bool, // True for right-hand twist, false for left-hand
    pub shooting_angle: f64,  // uphill/downhill angle in radians
    // MBA-717: real bullet geometry so spin-drift / Magnus / stability on this fast/MC
    // path use the actual bullet instead of hardcoded .308 / 1.24in / 10-twist placeholders.
    pub bullet_diameter: f64,                              // meters
    pub bullet_length: f64, // meters (0.0 -> derivatives falls back to the 4.5-caliber heuristic)
    pub twist_rate: f64,    // inches per turn
    pub custom_drag_table: Option<crate::drag::DragTable>, // Custom Drag Model (CDM) data
    /// MBA-1356: whole-curve scale for the custom deck (1.0 = neutral). Threaded so the
    /// binding entry point (fast_integrate_with_segments) cannot silently drop it.
    pub cd_scale: f64,
    /// Whether `custom_drag_table` is a reference curve that the BCs were measured against,
    /// flown as `Cd / BC` with the BC the kernel resolves (a ladder's band BC, else the scalar
    /// `bc`), rather than the projectile's own Cd over sectional density. Mirrors
    /// `BallisticInputs::custom_drag_reference` and is threaded beside `cd_scale` for the same
    /// reason: `build_inputs` copies it onto the derivative inputs, and a flag dropped beside
    /// `cd_scale = 1.0` would fly a reference curve as the projectile's own Cd. `false` for a
    /// caller with no table, or with the projectile's own curve.
    pub custom_drag_reference: bool,
    pub bc_segments: Option<Vec<(f64, f64)>>, // Mach-based BC segments: (mach, bc)
    pub use_bc_segments: bool,                // Whether to use BC segment interpolation
    /// MBA-1602: the caller's velocity-keyed BC ladder (fps bands), flown when
    /// `use_bc_segments` is on and the ladder is non-empty. `None` or an empty ladder is no
    /// ladder: with `use_bc_segments` on and no Mach `bc_segments`, the kernel keeps its
    /// historical opt-in auto-estimate (MBA-955). Threaded so the binding entry point
    /// (fast_integrate_with_segments) cannot silently drop it.
    pub bc_segments_data: Option<Vec<crate::BCSegmentData>>,
    /// MBA-954: altitude (m, relative to launch) below which integration stops. -1000.0 is the
    /// historical default — effectively "no early ground impact" for normal flat-fire shots.
    pub ground_threshold: f64,
    /// MBA-1137: optional downrange-segmented atmosphere. When `Some`, `compute_derivatives`
    /// swaps the standard-mode base T/P/H for the zone selected by downrange distance before the
    /// altitude lapse. `None` (default) is byte-identical to pre-feature behavior.
    pub atmo_sock: Option<crate::atmosphere::AtmoSock>,
}

/// Build the loop-invariant BallisticInputs for the derivatives function ONCE per integration,
/// instead of rebuilding it (a "none".to_string() alloc plus bc_segments / custom_drag_table
/// clones) on every derivative evaluation (4x per RK4 step, 7x per RK45 step). The launch-speed
/// magnitude supplies the muzzle-set Magnus spin; every other field depends only on `params`, so
/// the struct is constant for the whole integration.
fn build_inputs(params: &TrajectoryParams, muzzle_velocity_mps: f64) -> BallisticInputs {
    let mut inputs = BallisticInputs {
        bc_value: params.bc,
        bc_type: params.drag_model,
        // This generic RK4/RK45 kernel is fed a raw `TrajectoryParams.bc` that its callers
        // (fast_trajectory::fast_integrate_with_segments) read directly from
        // `BallisticInputs.bc_value` WITHOUT going through `TrajectorySolver::new` — the
        // single normalization boundary for MBA-1365. Any ASM-reference conversion is
        // therefore the caller's responsibility before populating `TrajectoryParams`;
        // this constant-for-the-whole-integration struct always reports ICAO (a no-op)
        // since it takes no reference-standard input of its own.
        bc_reference_standard: crate::cli_api::BcReferenceStandard::Icao,
        bullet_mass: params.mass_kg, // kg
        muzzle_velocity: muzzle_velocity_mps,
        bullet_diameter: params.bullet_diameter, // MBA-717: real geometry, not placeholders
        bullet_length: params.bullet_length,
        twist_rate: params.twist_rate,
        is_twist_right: params.is_twist_right,
        enable_advanced_effects: params.enable_spin_drift
            || params.enable_magnus
            || params.enable_coriolis,
        enable_magnus: params.enable_magnus,
        enable_coriolis: params.enable_coriolis,
        altitude: params.atmos_params.0,
        temperature: params.atmos_params.1,
        pressure: params.atmos_params.2,
        humidity: params.atmos_params.3,
        tipoff_yaw: 0.0,
        cd_delta2: 7.5,
        target_distance: 1000.0, // default
        muzzle_angle: 0.0,
        wind_speed: if !params.wind_segments.is_empty() {
            params.wind_segments[0].speed_kmh * 0.2777778 // km/h -> m/s
        } else {
            0.0
        },
        wind_angle: if !params.wind_segments.is_empty() {
            params.wind_segments[0].angle_deg.to_radians() // degrees -> radians
        } else {
            0.0
        },
        latitude: None,
        shooting_angle: params.shooting_angle,
        cant_angle: 0.0,
        azimuth_angle: 0.0,
        shot_azimuth: 0.0, // this fast path doesn't plumb latitude/bearing (no directional Coriolis here)
        use_powder_sensitivity: false,
        powder_temp_sensitivity: 0.0,
        powder_temp: 59.0,
        powder_temp_curve: None,
        powder_curve_temp_c: None,
        tipoff_decay_distance: 0.0,
        ground_threshold: params.ground_threshold, // MBA-954: honor the configured ground plane
        bc_segments: params.bc_segments.clone(),
        caliber_inches: params.bullet_diameter / 0.0254, // MBA-717: from real diameter
        weight_grains: params.mass_kg / crate::constants::GRAINS_TO_KG,
        use_bc_segments: params.use_bc_segments,
        bullet_id: None,
        // MBA-1602: fly the caller's ladder. An empty ladder is no ladder, so it collapses to
        // None here and the MBA-955 pre-population below behaves exactly as before.
        bc_segments_data: params
            .bc_segments_data
            .clone()
            .filter(|segments| !segments.is_empty()),
        use_enhanced_spin_drift: params.enable_spin_drift,
        use_form_factor: false,
        manufacturer: None,
        bullet_model: None,
        enable_wind_shear: false,
        wind_shear_model: "none".to_string(),
        use_cluster_bc: false,
        bullet_cluster: None,
        custom_drag_table: params.custom_drag_table.clone(),
        cd_scale: params.cd_scale,
        custom_drag_reference: params.custom_drag_reference,
        bc_type_str: None,
        enable_pitch_damping: false,
        enable_precession_nutation: false,
        // MBA-959/MBA-1183: aerodynamic jump stays OFF inside this low-level raw-state integrator.
        // The high-level fast wrappers form Sg from their complete BallisticInputs and rotate the
        // prebuilt initial velocity before entering their integration loops; enabling it again
        // here would double-apply the launch perturbation. Direct low-level callers likewise own
        // any desired launch-state rotation. (Real geometry is still carried for spin/Magnus.)
        enable_aerodynamic_jump: false,
        use_rk4: true,
        use_adaptive_rk45: false,
        enable_trajectory_sampling: false,
        sample_interval: 10.0,
        sight_height: 0.0,
        // MBA-1396: like cant_angle above, the fast kernel does not plumb sight-mount
        // geometry — callers own the initial lateral state.
        sight_offset_lateral_m: 0.0,
        muzzle_height: 0.0,
        target_height: 0.0,
        // MBA-1359: zero-POI offsets bias a ZERO SOLVE; this low-level kernel never solves
        // a zero (callers hand it a finished launch state), so they are structurally inert
        // here — same reasoning as cant_angle/shot_azimuth above.
        zero_poi_vertical_m: 0.0,
        zero_poi_horizontal_m: 0.0,
        // MBA-1403: an output-mode toggle for the trajectory sampler; this kernel never
        // runs trajectory sampling (enable_trajectory_sampling is false above), so the
        // default LOS reference is structurally inert here.
        drops_reference: crate::cli_api::DropsReference::Los,
    };

    // MBA-955: pre-populate velocity-BC segments ONCE here, instead of get_bc_for_velocity
    // rebuilding them (a model String + a segment Vec) on every derivative evaluation (4-7x per
    // step). Gated to EXACTLY the case where the per-step path would estimate: use_bc_segments on,
    // no explicit velocity segments, and no Mach-based bc_segments (those take a different,
    // unchanged path). bc_used there == params.bc == inputs.bc_value, so the estimated segments are
    // identical and the per-step fast-path lookup returns the same BC -> byte-identical output.
    if inputs.use_bc_segments && inputs.bc_segments_data.is_none() && inputs.bc_segments.is_none() {
        inputs.bc_segments_data =
            crate::derivatives::estimate_bc_segments_for(&inputs, inputs.bc_value);
    }
    inputs
}

/// Convert state to Vector6 and call compute_derivatives
fn compute_derivatives_vec_with_gravity(
    state: &Vector6<f64>,
    t: f64,
    params: &TrajectoryParams,
    inputs: &BallisticInputs,
    gravity: Gravity,
) -> Vector6<f64> {
    let pos = Vector3::new(state[0], state[1], state[2]);
    let vel = Vector3::new(state[3], state[4], state[5]);

    // Calculate wind at current position with shear support. Boundary-layer shear reads the
    // height above the ground: the shot-frame y, plus 0 when y is measured from the ground the
    // caller's `muzzle_height` was measured from, or the assumed bore height when y = 0 is the
    // launch point (`inputs.muzzle_height` 0, which every direct TrajectoryParams caller gets;
    // see integrate_trajectory_with_gravity). The ground is the one the ground_threshold check
    // uses, parallel to the line of sight, so the shooting angle does not enter. Site elevation
    // (`shooter_altitude_m`) is MSL and does not enter shear.
    let wind_vector = if !params.wind_segments.is_empty() {
        if params.enable_wind_shear && params.wind_shear_model != "none" {
            crate::wind_shear::apply_boundary_layer_shear_agl(
                wind_vector_for_range(pos.x, &params.wind_segments),
                crate::wind_shear::shear_height_agl_m(pos.y, inputs.muzzle_height),
                crate::wind_shear::boundary_layer_model_from_name(&params.wind_shear_model),
            )
        } else {
            wind_vector_for_range(pos.x, &params.wind_segments)
        }
    } else {
        Vector3::zeros()
    };

    // Call compute_derivatives - returns [f64; 6] directly. `inputs` is built once per
    // integration by build_inputs() and threaded in, instead of rebuilt every call.
    let deriv_result = compute_derivatives_with_gravity(
        pos,
        vel,
        inputs,
        wind_vector,
        params.atmos_params,
        params.bc,
        params.omega_vector,
        t,
        params.atmo_sock.as_ref(),
        gravity,
    );

    Vector6::new(
        deriv_result[0],
        deriv_result[1],
        deriv_result[2],
        deriv_result[3],
        deriv_result[4],
        deriv_result[5],
    )
}

/// Linearly localize a target crossing within an accepted forward integration step.
///
/// Callers provide a bracket with `start[0] <= target_x <= end[0]` and increasing downrange X.
/// The same crossing fraction is applied to time and every phase-space component; X is then set
/// exactly to the public target value to remove interpolation roundoff.
fn interpolate_target_crossing(
    start_time: f64,
    start: &Vector6<f64>,
    step_dt: f64,
    end: &Vector6<f64>,
    target_x: f64,
) -> (f64, Vector6<f64>) {
    debug_assert!(start[0] <= target_x && target_x <= end[0] && end[0] > start[0]);

    let alpha = (target_x - start[0]) / (end[0] - start[0]);
    let crossing_time = start_time + alpha * step_dt;
    let mut crossing_state = start + alpha * (end - start);
    crossing_state[0] = target_x;

    (crossing_time, crossing_state)
}

/// Checked sibling of [`integrate_trajectory`] (MBA-1338): rejects malformed wind
/// segments with a typed [`WindSegmentError`] **before** sorting, vector precomputation,
/// or producing any trajectory points, so a caller can never receive a poisoned
/// trajectory from non-finite segment fields. The error's `index` refers to the
/// caller's own segment ordering.
pub fn try_integrate_trajectory(
    initial_state: [f64; 6],
    t_span: (f64, f64),
    params: TrajectoryParams,
    method: &str,
    tolerance: f64,
    max_step: f64,
) -> Result<Vec<(f64, Vector6<f64>)>, WindSegmentError> {
    crate::wind::validate_wind_segments(&params.wind_segments)?;
    Ok(integrate_trajectory(
        initial_state,
        t_span,
        params,
        method,
        tolerance,
        max_step,
    ))
}

/// Main trajectory integration function
///
/// **Legacy/unchecked entry point** (MBA-1338): wind segments are sorted and consumed
/// without validation, so non-finite fields silently poison the returned trajectory.
/// Prefer [`try_integrate_trajectory`], which rejects malformed segments with a typed
/// error before any integration work.
pub fn integrate_trajectory(
    initial_state: [f64; 6],
    t_span: (f64, f64),
    params: TrajectoryParams,
    method: &str,
    tolerance: f64,
    max_step: f64,
) -> Vec<(f64, Vector6<f64>)> {
    integrate_trajectory_with_gravity(
        initial_state,
        t_span,
        params,
        method,
        tolerance,
        max_step,
        Gravity::EARTH,
    )
}

/// [`integrate_trajectory`] with an explicit gravity. `initial_state`'s y = 0 is the launch
/// point, as for every direct [`TrajectoryParams`] caller; see
/// [`integrate_trajectory_from_muzzle_height`] for a caller whose y is measured from the
/// ground.
#[allow(clippy::too_many_arguments)]
pub(crate) fn integrate_trajectory_with_gravity(
    initial_state: [f64; 6],
    t_span: (f64, f64),
    params: TrajectoryParams,
    method: &str,
    tolerance: f64,
    max_step: f64,
    gravity: Gravity,
) -> Vec<(f64, Vector6<f64>)> {
    integrate_trajectory_from_muzzle_height(
        initial_state,
        t_span,
        params,
        method,
        tolerance,
        max_step,
        gravity,
        0.0,
    )
}

/// [`integrate_trajectory_with_gravity`] for a launch state whose y is measured from the level
/// `muzzle_height_m` is measured from, as `BallisticInputs::muzzle_height` states it: the
/// ground when it is nonzero, the launch point (the bore) when it is 0. Only boundary-layer
/// shear reads it, to measure the height above the ground; this kernel used to add an assumed
/// 1.5 m bore height to every y, which counted a real bore height twice.
#[allow(clippy::too_many_arguments)]
pub(crate) fn integrate_trajectory_from_muzzle_height(
    initial_state: [f64; 6],
    t_span: (f64, f64),
    mut params: TrajectoryParams,
    method: &str,
    tolerance: f64,
    max_step: f64,
    gravity: Gravity,
    muzzle_height_m: f64,
) -> Vec<(f64, Vector6<f64>)> {
    // Normalize once before build_inputs reads the first zone and before any RK stage performs a
    // first-match lookup. Callers may supply zones in any order.
    crate::wind::sort_wind_segments_by_distance(&mut params.wind_segments);

    let mut state = Vector6::new(
        initial_state[0],
        initial_state[1],
        initial_state[2],
        initial_state[3],
        initial_state[4],
        initial_state[5],
    );

    let mut t = t_span.0;
    let t_end = t_span.1;
    let mut dt = (t_end - t) / 1000.0; // Initial step size

    let mut trajectory = Vec::with_capacity(10000);
    trajectory.push((t, state));
    if state[0] >= params.target_distance_m {
        return trajectory;
    }

    // Build the (loop-invariant) derivative inputs once for the whole integration, instead of
    // rebuilding the struct on every derivative evaluation.
    let muzzle_velocity_mps =
        Vector3::new(initial_state[3], initial_state[4], initial_state[5]).norm();
    let mut inputs = build_inputs(&params, muzzle_velocity_mps);
    // The frame datum is the caller's, not a TrajectoryParams field (adding one would break
    // every struct-literal constructor). The derivative model ignores muzzle_height; only the
    // shear height in compute_derivatives_vec_with_gravity reads it.
    inputs.muzzle_height = muzzle_height_m;

    match method {
        "RK4" => {
            // Fixed step RK4 with target detection
            dt = dt.min(max_step).min(0.001); // Use smaller steps for accuracy

            while t < t_end {
                if t + dt > t_end {
                    dt = t_end - t;
                }

                let new_state = rk4_step_with_gravity(&state, t, dt, &params, &inputs, gravity);

                // Check if we're about to pass the target (X is downrange, McCoy)
                if state[0] < params.target_distance_m && new_state[0] >= params.target_distance_m {
                    trajectory.push(interpolate_target_crossing(
                        t,
                        &state,
                        dt,
                        &new_state,
                        params.target_distance_m,
                    ));
                    break; // Stop at target
                }

                state = new_state;
                t += dt;
                trajectory.push((t, state));

                // Check if we've reached or passed the target
                if state[0] >= params.target_distance_m {
                    break;
                }

                // Check if bullet hit ground (MBA-954: honor the configured ground plane,
                // not a hardcoded -1000.0)
                if state[1] < params.ground_threshold {
                    break;
                }
            }
        }
        _ => {
            // Adaptive RK45 with better sampling
            let mut last_save_x = 0.0; // X is downrange (McCoy)
            let save_interval_m = params.target_distance_m / 50.0; // Save ~50 points minimum
            let tolerance = if tolerance.is_finite() && tolerance > 0.0 {
                tolerance
            } else {
                eprintln!(
                    "WARNING: RK45 tolerance must be finite and positive; using {RK45_DEFAULT_TOLERANCE}"
                );
                RK45_DEFAULT_TOLERANCE
            };

            // OPTIMIZATION: Adjust max step size when wind shear is enabled
            // This improves numerical stability at long ranges
            let effective_max_step =
                if params.enable_wind_shear && params.wind_shear_model != "none" {
                    // Use smaller steps for wind shear, but not TOO small
                    if params.target_distance_m > 800.0 {
                        0.01 // Smaller steps for long range with shear (10ms)
                    } else {
                        0.02 // Normal steps for medium range with shear (20ms)
                    }
                } else {
                    max_step // Use provided max_step when no wind shear
                };
            if !effective_max_step.is_finite() || effective_max_step <= 0.0 {
                eprintln!("WARNING: RK45 max_step must be finite and positive");
                return trajectory;
            }
            let min_step = RK45_MIN_STEP.min(effective_max_step);

            // Set initial step size - ensure it's reasonable
            dt = dt.min(effective_max_step).max(min_step);

            // Safety check: maximum iterations to prevent infinite loops
            let max_iterations = 100000; // Should be more than enough for any realistic trajectory
            let mut iteration_count = 0;

            while t < t_end && iteration_count < max_iterations {
                // Limit time step for better resolution
                if t + dt > t_end {
                    dt = t_end - t;
                }

                let control = Rk45Control {
                    tolerance,
                    min_step,
                    max_step: effective_max_step,
                    max_trials: max_iterations - iteration_count,
                };
                let accepted = match adaptive_rk45_step_with_gravity(
                    &state, t, dt, &params, &inputs, control, gravity,
                ) {
                    Ok(accepted) => accepted,
                    Err(trials) => {
                        iteration_count += trials;
                        if iteration_count < max_iterations {
                            eprintln!("WARNING: RK45 minimum-step trial was non-finite");
                        }
                        break;
                    }
                };
                iteration_count += accepted.trials;
                debug_assert!(accepted.error <= tolerance || accepted.used_dt <= min_step);

                // Target detection only examines an accepted candidate.
                if state[0] < params.target_distance_m
                    && accepted.state[0] >= params.target_distance_m
                {
                    trajectory.push(interpolate_target_crossing(
                        t,
                        &state,
                        accepted.used_dt,
                        &accepted.state,
                        params.target_distance_m,
                    ));
                    break;
                }

                // Update state and time using the interval that actually passed acceptance.
                state = accepted.state;
                t += accepted.used_dt;

                // Save trajectory point if we've moved enough distance
                if state[0] - last_save_x >= save_interval_m || state[0] >= params.target_distance_m
                {
                    // X is downrange
                    trajectory.push((t, state));
                    last_save_x = state[0];
                }

                // Limit the proposal for the next trial; this does not change the time just used.
                dt = accepted.next_dt;

                // Stop if we've reached the target
                if state[0] >= params.target_distance_m {
                    break;
                }

                // Check if bullet hit ground (MBA-954: honor the configured ground plane,
                // not a hardcoded -1000.0)
                if state[1] < params.ground_threshold {
                    break;
                }
            }

            // Warn if we hit the iteration limit
            if iteration_count >= max_iterations
                && t < t_end
                && state[0] < params.target_distance_m
                && state[1] >= params.ground_threshold
            {
                eprintln!(
                    "WARNING: Trajectory integration hit maximum iteration limit ({} iterations)",
                    max_iterations
                );
                eprintln!("  Final time: {}, Target time: {}", t, t_end);
                eprintln!(
                    "  Final position: downrange(x)={}, Target: {}m",
                    state[0], params.target_distance_m
                );
            }
        }
    }

    trajectory
}

/// Checked sibling of [`solve_trajectory_rust`] (MBA-1338): rejects malformed wind
/// segments with a typed [`WindSegmentError`] before any integration work or trajectory
/// points are produced. Bindings should migrate to this entry point so malformed
/// segments surface as a structured error instead of a silently poisoned trajectory.
#[allow(clippy::too_many_arguments)] // Mirrors the binding-compatibility signature below.
pub fn try_solve_trajectory_rust(
    initial_state: [f64; 6],
    t_span: (f64, f64),
    mass_kg: f64,
    bc: f64,
    drag_model: DragModel,
    wind_segments: Vec<WindSegment>,
    atmos_params: (f64, f64, f64, f64),
    omega_vector: Option<Vec<f64>>,
    enable_spin_drift: bool,
    enable_magnus: bool,
    enable_coriolis: bool,
    method: String,
    tolerance: f64,
    max_step: f64,
    target_distance_m: f64,
) -> Result<Vec<HashMap<String, f64>>, WindSegmentError> {
    crate::wind::validate_wind_segments(&wind_segments)?;
    Ok(solve_trajectory_rust(
        initial_state,
        t_span,
        mass_kg,
        bc,
        drag_model,
        wind_segments,
        atmos_params,
        omega_vector,
        enable_spin_drift,
        enable_magnus,
        enable_coriolis,
        method,
        tolerance,
        max_step,
        target_distance_m,
    ))
}

/// Python-exposed function for complete trajectory integration
///
/// **Legacy/unchecked entry point** (MBA-1338): consumes wind segments without
/// validation. Prefer [`try_solve_trajectory_rust`].
#[allow(clippy::too_many_arguments)] // Binding compatibility API; grouping would be breaking.
pub fn solve_trajectory_rust(
    initial_state: [f64; 6],
    t_span: (f64, f64),
    mass_kg: f64,
    bc: f64,
    drag_model: DragModel,
    wind_segments: Vec<WindSegment>,
    atmos_params: (f64, f64, f64, f64),
    omega_vector: Option<Vec<f64>>,
    enable_spin_drift: bool,
    enable_magnus: bool,
    enable_coriolis: bool,
    method: String,
    tolerance: f64,
    max_step: f64,
    target_distance_m: f64,
) -> Vec<HashMap<String, f64>> {
    let omega_vec = omega_vector.map(|v| Vector3::new(v[0], v[1], v[2]));

    let params = TrajectoryParams {
        mass_kg,
        bc,
        drag_model,
        wind_segments,
        atmos_params,
        omega_vector: omega_vec,
        enable_spin_drift,
        enable_magnus,
        enable_coriolis,
        target_distance_m,
        enable_wind_shear: false, // Default for test function
        wind_shear_model: "none".to_string(),
        shooter_altitude_m: 0.0,
        is_twist_right: true, // Default for test function
        shooting_angle: 0.0,  // This legacy entry takes no inclined-fire arg; flat fire only
        // This legacy entry takes no geometry args; keep the historical placeholders so its
        // behavior is unchanged (callers needing real geometry use fast_integrate_with_segments).
        bullet_diameter: 0.0078232,
        bullet_length: 0.031496,
        twist_rate: 10.0,
        custom_drag_table: None, // No CDM for test function
        cd_scale: 1.0,
        custom_drag_reference: false,
        bc_segments: None, // No BC segments for legacy function
        use_bc_segments: false,
        bc_segments_data: None,
        ground_threshold: -1000.0, // MBA-954: preserve the historical default
        atmo_sock: None,           // MBA-1137: legacy entry has no downrange atmosphere
    };

    let trajectory =
        integrate_trajectory(initial_state, t_span, params, &method, tolerance, max_step);

    // Convert to Python-friendly format
    trajectory
        .into_iter()
        .map(|(t, state)| {
            let mut point = HashMap::new();
            point.insert("t".to_string(), t);
            point.insert("x".to_string(), state[0]);
            point.insert("y".to_string(), state[1]);
            point.insert("z".to_string(), state[2]);
            point.insert("vx".to_string(), state[3]);
            point.insert("vy".to_string(), state[4]);
            point.insert("vz".to_string(), state[5]);
            point
        })
        .collect()
}

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

    #[test]
    fn uniform_gravity_integrators_match_a_ball_in_near_vacuum() {
        let initial = [0.0, 10.0, 0.0, 10.0, 2.0, 0.0];
        for method in ["RK4", "RK45"] {
            for magnitude in [0.0, 1.62, 9.80665, 20.0] {
                let mut params = create_test_params(10.0);
                params.mass_kg = 0.0459;
                params.bullet_diameter = 0.0427;
                params.bullet_length = 0.0427;
                params.atmos_params = (1e-12, 340.0, 0.0, 0.0);
                let result = integrate_trajectory_with_gravity(
                    initial,
                    (0.0, 2.0),
                    params,
                    method,
                    1e-8,
                    0.001,
                    Gravity::new(magnitude).unwrap(),
                );
                let (time, state) = result.last().unwrap();
                assert!((time - 1.0).abs() < 1e-8);
                assert!(
                    (state[1] - (10.0 + 2.0 * time - 0.5 * magnitude * time * time)).abs() < 1e-5
                );
                assert!((state[4] - (2.0 - magnitude * time)).abs() < 1e-8);
            }
        }
    }

    fn create_test_params(target_distance_m: f64) -> TrajectoryParams {
        TrajectoryParams {
            mass_kg: 0.01134, // 175 grains in kg
            bc: 0.442,
            bullet_diameter: 0.0078232, // .308 in
            bullet_length: 0.031496,    // 1.24 in
            twist_rate: 10.0,
            drag_model: DragModel::G7,
            wind_segments: vec![],
            atmos_params: (0.0, 15.0, 1013.25, 1.0),
            omega_vector: None,
            enable_spin_drift: false,
            enable_magnus: false,
            enable_coriolis: false,
            target_distance_m,
            enable_wind_shear: false,
            wind_shear_model: "none".to_string(),
            shooter_altitude_m: 0.0,
            is_twist_right: true,
            shooting_angle: 0.0,
            custom_drag_table: None,
            cd_scale: 1.0,
            custom_drag_reference: false,
            bc_segments: None,
            use_bc_segments: false,
            bc_segments_data: None,
            ground_threshold: -1000.0,
            atmo_sock: None,
        }
    }

    /// `build_inputs` is where the derivative kernel's inputs are rebuilt from the params, so a
    /// field it forgets is silently dropped (MBA-1356). Both custom-table fields arrive as set.
    #[test]
    fn build_inputs_carries_the_custom_table_fields() {
        for (cd_scale, reference) in [(1.0, false), (1.25, false), (1.0, true), (0.8, true)] {
            let params = TrajectoryParams {
                custom_drag_table: Some(crate::drag::DragTable::new(
                    vec![0.5, 3.0],
                    vec![0.2, 0.3],
                )),
                cd_scale,
                custom_drag_reference: reference,
                ..create_test_params(300.0)
            };
            let inputs = build_inputs(&params, 800.0);
            assert!(inputs.custom_drag_table.is_some());
            assert_eq!(inputs.cd_scale, cd_scale);
            assert_eq!(inputs.custom_drag_reference, reference);
        }
    }

    #[test]
    fn try_integrate_trajectory_rejects_malformed_segments_before_any_points() {
        // MBA-1338: a poisoned segment (index 1, caller order) must yield the typed error
        // and NO trajectory points — validation precedes sorting and integration.
        let mut params = create_test_params(300.0);
        params.wind_segments = vec![
            WindSegment::new(10.0, 90.0, 200.0),
            WindSegment::new(10.0, 90.0, f64::NAN),
        ];
        let err = try_integrate_trajectory(
            [0.0, 0.0, 0.0, 800.0, 0.0, 0.0],
            (0.0, 2.0),
            params,
            "RK4",
            1e-6,
            0.001,
        )
        .unwrap_err();
        assert_eq!(err.index, 1);
        assert_eq!(err.field, crate::wind::WindSegmentField::UntilM);
        assert_eq!(
            err.to_string(),
            "wind.segments[1].until_m must be finite and greater than zero"
        );
    }

    #[test]
    fn try_integrate_trajectory_matches_unchecked_on_valid_input() {
        let mk = || {
            let mut params = create_test_params(300.0);
            params.wind_segments = vec![WindSegment::new(16.0934, 90.0, 500.0)];
            params
        };
        let checked = try_integrate_trajectory(
            [0.0, 0.0, 0.0, 800.0, 0.0, 0.0],
            (0.0, 2.0),
            mk(),
            "RK4",
            1e-6,
            0.001,
        )
        .expect("valid segments must integrate");
        let unchecked = integrate_trajectory(
            [0.0, 0.0, 0.0, 800.0, 0.0, 0.0],
            (0.0, 2.0),
            mk(),
            "RK4",
            1e-6,
            0.001,
        );
        assert_eq!(checked.len(), unchecked.len());
        assert_eq!(checked.last().unwrap().1, unchecked.last().unwrap().1);
    }

    #[test]
    fn try_solve_trajectory_rust_rejects_malformed_segments() {
        let bad = vec![WindSegment::new(-5.0, 0.0, 100.0)];
        let err = try_solve_trajectory_rust(
            [0.0, 0.0, 0.0, 800.0, 0.0, 0.0],
            (0.0, 2.0),
            0.01134,
            0.442,
            DragModel::G7,
            bad,
            (0.0, 15.0, 1013.25, 1.0),
            None,
            false,
            false,
            false,
            "RK4".to_string(),
            1e-6,
            0.001,
            300.0,
        )
        .unwrap_err();
        assert_eq!(err.index, 0);
        assert_eq!(err.field, crate::wind::WindSegmentField::SpeedKmh);
        assert_eq!(
            err.to_string(),
            "wind.segments[0].speed_kmh must be finite and non-negative"
        );
    }

    #[test]
    fn try_solve_trajectory_rust_succeeds_on_valid_segments() {
        let points = try_solve_trajectory_rust(
            [0.0, 0.0, 0.0, 800.0, 0.0, 0.0],
            (0.0, 2.0),
            0.01134,
            0.442,
            DragModel::G7,
            vec![WindSegment::new(16.0934, 90.0, 500.0)],
            (0.0, 15.0, 1013.25, 1.0),
            None,
            false,
            false,
            false,
            "RK4".to_string(),
            1e-6,
            0.001,
            300.0,
        )
        .expect("valid segments must solve");
        assert!(!points.is_empty());
        assert!(points.last().unwrap()["x"] > 0.0);
    }

    #[test]
    fn derivative_inputs_preserve_initial_velocity_as_muzzle_speed() {
        let params = create_test_params(1_000.0);
        let launch_velocity: Vector3<f64> = Vector3::new(700.0, 30.0, -20.0);
        let inputs = build_inputs(&params, launch_velocity.norm());

        assert_eq!(
            inputs.muzzle_velocity.to_bits(),
            launch_velocity.norm().to_bits()
        );
    }

    #[test]
    fn integrated_magnus_retains_nonzero_launch_spin() {
        let initial_state = [0.0, 0.0, 0.0, 800.0, 0.0, 0.0];
        let baseline = integrate_trajectory(
            initial_state,
            (0.0, 0.1),
            create_test_params(1_000.0),
            "RK4",
            1e-6,
            0.001,
        );
        let mut magnus_params = create_test_params(1_000.0);
        magnus_params.enable_magnus = true;

        let trajectory =
            integrate_trajectory(initial_state, (0.0, 0.1), magnus_params, "RK4", 1e-6, 0.001);
        let baseline_y = baseline.last().expect("baseline trajectory is empty").1[1];
        let magnus_y = trajectory.last().expect("trajectory is empty").1[1];
        let vertical_delta = magnus_y - baseline_y;

        assert!(
            vertical_delta.is_finite() && vertical_delta < 0.0,
            "right-twist Magnus should retain nonzero launch spin and point down, got \
             delta_y={vertical_delta}"
        );
    }

    #[test]
    fn rk45_retries_rejected_wind_boundary_step() {
        let initial_state = [0.0, 0.0, 0.0, 800.0, 0.0, 0.0];
        let mut params = create_test_params(100.0);
        params.wind_segments = vec![
            WindSegment::new(0.0, 90.0, 4.0),
            WindSegment::new(1_000.0, 90.0, 10_000.0),
        ];

        let state = Vector6::from_row_slice(&initial_state);
        let launch_speed =
            Vector3::new(initial_state[3], initial_state[4], initial_state[5]).norm();
        let inputs = build_inputs(&params, launch_speed);
        let initial_dt = 0.01;
        let tolerance = 1e-6;
        let (rejected_state, suggested_dt, error) = rk45_step_with_gravity(
            &state,
            0.0,
            initial_dt,
            &params,
            &inputs,
            tolerance,
            Gravity::EARTH,
        );
        assert!(
            error > tolerance,
            "wind-boundary trial must exceed tolerance, got {error}"
        );
        assert!(suggested_dt < initial_dt);

        let accepted = adaptive_rk45_step_with_gravity(
            &state,
            0.0,
            initial_dt,
            &params,
            &inputs,
            Rk45Control {
                tolerance,
                min_step: RK45_MIN_STEP,
                max_step: initial_dt,
                max_trials: 100,
            },
            Gravity::EARTH,
        )
        .expect("a smaller finite trial should satisfy the tolerance");

        assert!(accepted.trials > 1, "oversized trial was not retried");
        assert!(accepted.used_dt < initial_dt);
        assert!(
            accepted.error <= tolerance || accepted.used_dt <= RK45_MIN_STEP,
            "accepted error {} exceeds tolerance at dt {}",
            accepted.error,
            accepted.used_dt
        );

        let (accepted_state, _, accepted_error) = rk45_step_with_gravity(
            &state,
            0.0,
            accepted.used_dt,
            &params,
            &inputs,
            tolerance,
            Gravity::EARTH,
        );
        assert_eq!(accepted.state, accepted_state);
        assert_eq!(accepted.error, accepted_error);
        assert_ne!(accepted.state, rejected_state);
        assert!((RK45_MIN_STEP..=initial_dt).contains(&accepted.next_dt));
    }

    #[test]
    fn integration_normalizes_wind_segments_by_distance() {
        let initial_state = [0.0, 0.0, 0.0, 800.0, 0.0, 0.0];
        let sorted_segments = vec![
            WindSegment::new(40.0, 270.0, 300.0),
            WindSegment::new(20.0, 90.0, 600.0),
        ];

        let mut sorted_params = create_test_params(100.0);
        sorted_params.wind_segments = sorted_segments.clone();
        let mut unsorted_params = create_test_params(100.0);
        unsorted_params.wind_segments = sorted_segments.into_iter().rev().collect();

        let sorted =
            integrate_trajectory(initial_state, (0.0, 1.0), sorted_params, "RK4", 1e-6, 0.001);
        let unsorted = integrate_trajectory(
            initial_state,
            (0.0, 1.0),
            unsorted_params,
            "RK4",
            1e-6,
            0.001,
        );

        assert_eq!(unsorted.len(), sorted.len());
        for (index, ((sorted_t, sorted_state), (unsorted_t, unsorted_state))) in
            sorted.iter().zip(&unsorted).enumerate()
        {
            assert_eq!(unsorted_t.to_bits(), sorted_t.to_bits());
            for component in 0..6 {
                assert_eq!(
                    unsorted_state[component].to_bits(),
                    sorted_state[component].to_bits(),
                    "wind segment order changed state component {component} at point {index}"
                );
            }
        }
    }

    #[test]
    fn rk4_target_crossing_interpolates_complete_state_and_time() {
        let initial_state = [0.0, 0.0, 0.0, 800.0, 5.0, 2.0];
        let target_distance_m = 100.0;
        let trajectory = integrate_trajectory(
            initial_state,
            (0.0, 1.0),
            create_test_params(target_distance_m),
            "RK4",
            1e-6,
            0.001,
        );

        let (previous_t, previous_state) = &trajectory[trajectory.len() - 2];
        let (terminal_t, terminal_state) = trajectory.last().expect("trajectory is empty");
        let reference_params = create_test_params(target_distance_m);
        let inputs = build_inputs(&reference_params, Vector3::new(800.0, 5.0, 2.0).norm());
        let full_step_dt = 0.001;
        let bracket_end = rk4_step_with_gravity(
            previous_state,
            *previous_t,
            full_step_dt,
            &reference_params,
            &inputs,
            Gravity::EARTH,
        );
        assert!(previous_state[0] < target_distance_m);
        assert!(bracket_end[0] >= target_distance_m);

        let alpha = (target_distance_m - previous_state[0]) / (bracket_end[0] - previous_state[0]);
        let expected_t = previous_t + alpha * full_step_dt;
        let mut expected_state = previous_state + alpha * (bracket_end - previous_state);
        expected_state[0] = target_distance_m;

        assert_eq!(terminal_t.to_bits(), expected_t.to_bits());
        for component in 0..6 {
            assert_eq!(
                terminal_state[component].to_bits(),
                expected_state[component].to_bits(),
                "terminal component {component} was not interpolated at the target crossing"
            );
        }
    }

    #[test]
    fn rk45_target_crossing_uses_the_accepted_state_and_time() {
        let initial_state = [0.0, 0.0, 0.0, 800.0, 5.0, 2.0];
        let initial = Vector6::from_row_slice(&initial_state);
        let target_distance_m = 0.5;
        let reference_params = create_test_params(target_distance_m);
        let inputs = build_inputs(&reference_params, Vector3::new(800.0, 5.0, 2.0).norm());
        let initial_dt = 0.001;
        let accepted = adaptive_rk45_step_with_gravity(
            &initial,
            0.0,
            initial_dt,
            &reference_params,
            &inputs,
            Rk45Control {
                tolerance: 1e-6,
                min_step: RK45_MIN_STEP,
                max_step: 0.01,
                max_trials: 100_000,
            },
            Gravity::EARTH,
        )
        .expect("first RK45 target bracket should be accepted");
        assert!(accepted.state[0] >= target_distance_m);
        let expected = interpolate_target_crossing(
            0.0,
            &initial,
            accepted.used_dt,
            &accepted.state,
            target_distance_m,
        );

        let trajectory = integrate_trajectory(
            initial_state,
            (0.0, 1.0),
            create_test_params(target_distance_m),
            "RK45",
            1e-6,
            0.01,
        );
        let actual = trajectory.last().expect("trajectory is empty");

        assert_eq!(actual.0.to_bits(), expected.0.to_bits());
        for component in 0..6 {
            assert_eq!(
                actual.1[component].to_bits(),
                expected.1[component].to_bits(),
                "RK45 terminal component {component} was not interpolated from its accepted step"
            );
        }
    }

    #[test]
    fn target_crossing_helper_interpolates_every_component() {
        let start = Vector6::new(90.0, 10.0, -4.0, 700.0, -20.0, 5.0);
        let end = Vector6::new(130.0, 6.0, 8.0, 660.0, -24.0, 9.0);
        let (time, state) = interpolate_target_crossing(2.0, &start, 0.5, &end, 100.0);

        assert_eq!(time.to_bits(), 2.125_f64.to_bits());
        for (index, expected) in [100.0_f64, 9.0, -1.0, 690.0, -21.0, 6.0]
            .into_iter()
            .enumerate()
        {
            assert_eq!(state[index].to_bits(), expected.to_bits());
        }
    }

    #[test]
    fn already_at_or_past_target_returns_initial_state_without_advancing() {
        let initial = [150.0, 12.0, -3.0, 700.0, -4.0, 5.0];

        for method in ["RK4", "RK45"] {
            for target in [150.0, 100.0] {
                let trajectory = integrate_trajectory(
                    initial,
                    (2.0, 3.0),
                    create_test_params(target),
                    method,
                    1e-6,
                    0.01,
                );

                assert_eq!(trajectory.len(), 1, "{method} advanced a terminal state");
                let (time, state) = &trajectory[0];
                assert_eq!(time.to_bits(), 2.0_f64.to_bits());
                for index in 0..6 {
                    assert_eq!(state[index].to_bits(), initial[index].to_bits());
                }
            }
        }
    }

    #[test]
    fn rk45_error_norm_scales_components_independently() {
        let state = Vector6::new(1.0e9, 0.0, 0.0, 800.0, 0.0, 0.0);
        let fifth_order = state;
        let mut fourth_order = state;
        fourth_order[4] = 1.0e-3;

        let error = rk45_error_norm(&state, &fifth_order, &fourth_order);
        let expected = 1.0e-3 / 6.0_f64.sqrt();

        assert!(
            (error - expected).abs() <= 1e-15,
            "large downrange position masked a velocity-component error: {error}"
        );
    }

    #[test]
    fn test_mba954_ground_threshold_honored() {
        // MBA-954: integrate_trajectory must honor the configured ground plane, not a hardcoded
        // -1000.0. A descending bullet with a shallow ground_threshold must terminate earlier
        // (fewer points) than one with the historical deep default.
        let initial_state = [0.0, 0.0, 0.0, 300.0, -30.0, 0.0]; // descending (vy = -30 m/s)

        let mut shallow = create_test_params(1_000_000.0); // huge target so range never terminates
        shallow.ground_threshold = -20.0; // stop ~20 m below launch
        let mut deep = create_test_params(1_000_000.0);
        deep.ground_threshold = -1000.0; // historical default

        let t_shallow =
            integrate_trajectory(initial_state, (0.0, 60.0), shallow, "RK4", 1e-6, 0.001);
        let t_deep = integrate_trajectory(initial_state, (0.0, 60.0), deep, "RK4", 1e-6, 0.001);

        assert!(
            t_shallow.len() < t_deep.len(),
            "shallow ground_threshold (-20) should terminate earlier than deep (-1000): \
             shallow={}, deep={}",
            t_shallow.len(),
            t_deep.len()
        );
    }

    #[test]
    fn test_integrate_trajectory_basic() {
        // Initial state [x,y,z,vx,vy,vz] (McCoy: X=downrange, Z=lateral)
        // x=0 (downrange start), vx=821.52 (downrange velocity)
        let initial_state = [0.0, -0.038, 0.0, 821.52, 48.61, 0.0];

        let params = TrajectoryParams {
            mass_kg: 0.01134, // 175 grains in kg
            bc: 0.442,
            bullet_diameter: 0.0078232, // .308 in
            bullet_length: 0.031496,    // 1.24 in
            twist_rate: 10.0,
            drag_model: DragModel::G7,
            wind_segments: vec![WindSegment::new(0.0, 90.0, 914.4)],
            atmos_params: (0.0, 15.0, 1013.25, 1.0),
            omega_vector: None,
            enable_spin_drift: false,
            enable_magnus: false,
            enable_coriolis: false,
            target_distance_m: 914.4, // 1000 yards in meters
            enable_wind_shear: false,
            wind_shear_model: "none".to_string(),
            shooter_altitude_m: 0.0,
            is_twist_right: true,
            shooting_angle: 0.0,
            custom_drag_table: None,
            cd_scale: 1.0,
            custom_drag_reference: false,
            bc_segments: None,
            use_bc_segments: false,
            bc_segments_data: None,
            ground_threshold: -1000.0,
            atmo_sock: None,
        };

        println!("Running integrate_trajectory test...");
        println!("Initial state: {:?}", initial_state);
        println!("Target distance: {} m", params.target_distance_m);

        let trajectory =
            integrate_trajectory(initial_state, (0.0, 10.0), params, "RK45", 1e-6, 0.01);

        println!("Trajectory has {} points", trajectory.len());

        // Should have more than just initial point
        assert!(
            trajectory.len() > 1,
            "Trajectory should have more than 1 point, but has {}",
            trajectory.len()
        );

        // Check that we actually moved downrange
        if let Some((_, final_state)) = trajectory.last() {
            println!("Final state: downrange(x)={}", final_state[0]);
            assert!(
                final_state[0] > 0.0,
                "Final x should be positive (bullet moved downrange)"
            );
            assert!(
                final_state[0] >= 900.0,
                "Final x should be near target distance"
            );
            assert!(
                final_state[3] < 0.9 * initial_state[3],
                "standard-atmosphere drag should reduce downrange velocity"
            );
        }
    }

    #[test]
    fn test_rk4_vs_rk45_consistency() {
        // Both methods should give similar results for the same trajectory
        let initial_state = [0.0, 0.0, 0.0, 800.0, 30.0, 0.0]; // McCoy: vx=downrange
        let target_distance = 500.0;

        let params_rk4 = create_test_params(target_distance);
        let params_rk45 = create_test_params(target_distance);

        let trajectory_rk4 =
            integrate_trajectory(initial_state, (0.0, 5.0), params_rk4, "RK4", 1e-6, 0.001);
        let trajectory_rk45 =
            integrate_trajectory(initial_state, (0.0, 5.0), params_rk45, "RK45", 1e-6, 0.01);

        // Both should reach target
        assert!(!trajectory_rk4.is_empty());
        assert!(!trajectory_rk45.is_empty());

        let (time_rk4, final_rk4) = trajectory_rk4.last().unwrap();
        let (time_rk45, final_rk45) = trajectory_rk45.last().unwrap();

        // Compare quantities that are not forced equal by target-distance clamping.
        assert!(
            (time_rk4 - time_rk45).abs() < 1e-4,
            "RK4/RK45 time of flight diverged: {time_rk4} vs {time_rk45}"
        );
        assert!((final_rk4[1] - final_rk45[1]).abs() < 1e-3);
        assert!((final_rk4[3] - final_rk45[3]).abs() < 1e-2);
        assert!(final_rk45[3] < 0.9 * initial_state[3]);
    }

    #[test]
    fn test_ground_impact_detection() {
        // Trajectory with steep downward angle should hit ground
        let initial_state = [0.0, 100.0, 0.0, 300.0, -50.0, 0.0]; // McCoy: vx=downrange // Steep descent

        let mut params = create_test_params(10000.0); // Far target
        params.target_distance_m = 10000.0;
        let ground_threshold = 0.0;
        params.ground_threshold = ground_threshold;

        let trajectory =
            integrate_trajectory(initial_state, (0.0, 20.0), params, "RK4", 1e-6, 0.01);

        // Should stop before reaching target due to ground impact
        let (_, final_state) = trajectory.last().unwrap();

        // y should have crossed the configured ground threshold.
        assert!(
            final_state[1] <= ground_threshold,
            "Should hit ground, but y={}",
            final_state[1]
        );
        assert!(
            final_state[0] < 10000.0,
            "Should not reach target, but z={}",
            final_state[0]
        );
    }

    #[test]
    fn test_target_distance_reached() {
        let initial_state = [0.0, 0.0, 0.0, 800.0, 20.0, 0.0]; // McCoy: vx=downrange
        let target_distance = 300.0;

        let params = create_test_params(target_distance);

        let trajectory =
            integrate_trajectory(initial_state, (0.0, 5.0), params, "RK45", 1e-6, 0.01);

        let (_, final_state) = trajectory.last().unwrap();

        // Should stop at or very near target distance
        assert!(
            (final_state[0] - target_distance).abs() < 1.0,
            "Should reach target at {}m, but stopped at {}m",
            target_distance,
            final_state[0]
        );
    }

    #[test]
    fn test_wind_affects_trajectory() {
        // Test that wind segments are properly stored and passed through
        // The actual wind effect depends on the derivatives computation which
        // uses the wind vector in the drag calculation
        let initial_state = [0.0, 0.0, 0.0, 800.0, 30.0, 0.0]; // McCoy: vx=downrange
        let target_distance = 500.0;

        // No wind
        let params_no_wind = create_test_params(target_distance);

        // Strong headwind (0 degrees = headwind)
        let mut params_headwind = create_test_params(target_distance);
        params_headwind.wind_segments = vec![WindSegment::new(72.0, 0.0, 500.0)]; // 72 km/h = 20 m/s headwind

        let trajectory_no_wind = integrate_trajectory(
            initial_state,
            (0.0, 5.0),
            params_no_wind,
            "RK45",
            1e-6,
            0.01,
        );
        let trajectory_headwind = integrate_trajectory(
            initial_state,
            (0.0, 5.0),
            params_headwind,
            "RK45",
            1e-6,
            0.01,
        );

        // Both trajectories should complete
        assert!(
            !trajectory_no_wind.is_empty(),
            "No-wind trajectory should complete"
        );
        assert!(
            !trajectory_headwind.is_empty(),
            "Headwind trajectory should complete"
        );

        let (time_no_wind, final_no_wind) = trajectory_no_wind.last().unwrap();
        let (time_headwind, final_headwind) = trajectory_headwind.last().unwrap();

        // Headwind should slow the bullet, resulting in longer flight time
        // or different drop at same distance
        let drop_no_wind = final_no_wind[1];
        let drop_headwind = final_headwind[1];

        println!("No wind: time={}, drop={}", time_no_wind, drop_no_wind);
        println!("Headwind: time={}, drop={}", time_headwind, drop_headwind);

        assert!(
            *time_headwind > *time_no_wind + 0.001,
            "headwind should increase time of flight: no-wind={time_no_wind}, headwind={time_headwind}"
        );
        assert!(
            final_headwind[3] < final_no_wind[3] - 1.0,
            "headwind should reduce terminal downrange velocity"
        );

        // Both should reach approximately the target distance
        assert!(
            (final_no_wind[0] - target_distance).abs() < 10.0,
            "No-wind should reach target"
        );
        assert!(
            (final_headwind[0] - target_distance).abs() < 10.0,
            "Headwind should reach target"
        );
    }

    #[test]
    fn test_solve_trajectory_rust_output_format() {
        let initial_state = [0.0, 0.0, 0.0, 800.0, 30.0, 0.0]; // McCoy: vx=downrange

        let result = solve_trajectory_rust(
            initial_state,
            (0.0, 2.0),
            0.01134,       // mass_kg
            0.442,         // bc
            DragModel::G7, // drag_model
            vec![],        // wind_segments
            // Standard atmosphere: altitude m, temperature C, pressure hPa, density ratio.
            (0.0, 15.0, 1013.25, 1.0),
            None,               // omega_vector
            false,              // enable_spin_drift
            false,              // enable_magnus
            false,              // enable_coriolis
            "RK45".to_string(), // method
            1e-6,               // tolerance
            0.01,               // max_step
            500.0,              // target_distance_m
        );

        // Should return Vec of HashMaps with expected keys
        assert!(!result.is_empty());

        let first_point = &result[0];
        assert!(first_point.contains_key("t"));
        assert!(first_point.contains_key("x"));
        assert!(first_point.contains_key("y"));
        assert!(first_point.contains_key("z"));
        assert!(first_point.contains_key("vx"));
        assert!(first_point.contains_key("vy"));
        assert!(first_point.contains_key("vz"));

        let final_point = result.last().unwrap();
        assert!(
            final_point["vx"] < 0.9 * initial_state[3],
            "standard-atmosphere wrapper fixture should exercise drag"
        );
    }

    #[test]
    fn test_left_vs_right_twist() {
        let initial_state = [0.0, 0.0, 0.0, 800.0, 30.0, 0.0]; // McCoy: vx=downrange
        let target_distance = 500.0;

        let mut params_right = create_test_params(target_distance);
        params_right.is_twist_right = true;
        params_right.enable_spin_drift = true;

        let mut params_left = create_test_params(target_distance);
        params_left.is_twist_right = false;
        params_left.enable_spin_drift = true;

        let trajectory_right =
            integrate_trajectory(initial_state, (0.0, 5.0), params_right, "RK45", 1e-6, 0.01);
        let trajectory_left =
            integrate_trajectory(initial_state, (0.0, 5.0), params_left, "RK45", 1e-6, 0.01);

        // Both should complete
        assert!(!trajectory_right.is_empty());
        assert!(!trajectory_left.is_empty());

        // Right and left twist should produce valid trajectories
        let (_, final_right) = trajectory_right.last().unwrap();
        let (_, final_left) = trajectory_left.last().unwrap();

        // Both should reach approximately the same downrange distance
        assert!((final_right[2] - final_left[2]).abs() < 10.0);
    }

    /// Asserts that the boundary-layer multiplier this kernel applies at `(x, y)` is
    /// `expected_ratio`: its derivatives with shear on equal its derivatives with shear off and
    /// the wind scaled by that ratio. `muzzle_height` is the level the caller's y is measured
    /// from, as `integrate_trajectory_with_gravity` sets it on the derivative inputs.
    fn assert_kernel_shear_ratio(
        model: &str,
        muzzle_height: f64,
        shooting_angle: f64,
        (x, y): (f64, f64),
        expected_ratio: f64,
    ) {
        let state = Vector6::new(x, y, 0.0, 700.0, 5.0, 0.0);
        let run = |enable_wind_shear: bool, wind_speed_kmh: f64| {
            let mut params = create_test_params(2_000.0);
            params.wind_segments = vec![WindSegment::new(wind_speed_kmh, 90.0, 2_000.0)];
            params.enable_wind_shear = enable_wind_shear;
            params.wind_shear_model = model.to_string();
            params.shooting_angle = shooting_angle;
            let mut inputs = build_inputs(&params, 700.0);
            inputs.muzzle_height = muzzle_height;
            compute_derivatives_vec_with_gravity(&state, 0.0, &params, &inputs, Gravity::EARTH)
        };
        let sheared = run(true, 36.0);
        let scaled = run(false, 36.0 * expected_ratio);
        for component in 3..6 {
            assert!(
                (sheared[component] - scaled[component]).abs() < 1e-12,
                "{model}, muzzle_height {muzzle_height} m, shooting angle {shooting_angle} rad, \
                 ({x}, {y}): the shear multiplier is not {expected_ratio} (component \
                 {component}: {} vs {})",
                sheared[component],
                scaled[component]
            );
        }
    }

    fn shear_ratio_rel_launch(height_rel_launch_m: f64, model: &str) -> f64 {
        crate::wind_shear::boundary_layer_speed_ratio(
            height_rel_launch_m,
            crate::wind_shear::boundary_layer_model_from_name(model),
        )
    }

    /// A direct `TrajectoryParams` caller's y = 0 is the launch point (the bore), so on a level
    /// shot the shear height is y plus the assumed 1.5 m bore height, as before. A caller that
    /// states a real `muzzle_height` (fast_integrate_with_segments forwards its inputs') has
    /// y = 0 on the ground instead, and the same height above the ground reads the same
    /// multiplier.
    #[test]
    fn kernel_shear_reads_one_multiplier_at_one_height_above_ground() {
        for model in ["logarithmic", "power_law"] {
            let expected = shear_ratio_rel_launch(58.5, model);
            assert!(
                expected > 1.2,
                "{model}: the fixture must be well above the floor"
            );
            assert_kernel_shear_ratio(model, 0.0, 0.0, (1_000.0, 58.5), expected);
            assert_kernel_shear_ratio(model, 1.5, 0.0, (1_000.0, 60.0), expected);
        }
    }

    /// On an inclined shot this kernel works in the shot frame like the other two, so its shear
    /// height is the shot-frame y plus the height of y = 0 above the ground, the ground running
    /// parallel to the line of sight as it does for the ground_threshold check, not the world
    /// height above the shooter's level, x·sin θ + y·cos θ.
    #[test]
    fn inclined_kernel_shear_measures_height_along_the_line_of_sight() {
        let (x, y) = (800.0, 60.0);
        for angle in [0.3_f64, -0.3] {
            let world_height = crate::atmosphere::shot_frame_altitude(0.0, x, y, angle);
            for model in ["logarithmic", "power_law"] {
                let bore_frame = shear_ratio_rel_launch(y, model);
                let ground_frame = shear_ratio_rel_launch(y - 1.5, model);
                assert!(
                    bore_frame > 1.25 && ground_frame > 1.25,
                    "{model}: fixture too low"
                );
                assert!(
                    (shear_ratio_rel_launch(world_height, model) - bore_frame).abs() > 0.1,
                    "{model}, angle {angle}: the fixture must tell the world height apart"
                );
                assert_kernel_shear_ratio(model, 0.0, angle, (x, y), bore_frame);
                assert_kernel_shear_ratio(model, 1.5, angle, (x, y), ground_frame);
            }
        }
    }
}