ballistics_engine/

derivatives.rs

use crate::atmosphere::{
    calculate_air_density_cimp, get_direct_atmosphere, get_local_atmosphere_humid, AtmoSock,
};
use crate::bc_estimation::{velocity_segment_bc, BCSegmentEstimator};
use crate::constants::*;
use crate::drag::get_drag_coefficient_full;
use crate::InternalBallisticInputs as BallisticInputs;
use nalgebra::Vector3;

// Magnus Effect Constants
//
// The Magnus effect causes spinning projectiles to deflect perpendicular to both
// their velocity vector and spin axis due to asymmetric pressure distribution.
// These constants define the Magnus moment coefficient (C_Lα) for different flight regimes.

/// Magnus coefficient for subsonic flow (M < 0.8)
///
/// Value: 0.030 (dimensionless coefficient)
/// Physical basis: Fully developed boundary layer circulation around spinning projectile
/// Regime: Subsonic flow where boundary layer remains attached
/// Source: McCoy's "Modern Exterior Ballistics", validated against wind tunnel data
const MAGNUS_COEFF_SUBSONIC: f64 = 0.030;

/// Magnus coefficient reduction factor for transonic regime (0.8 < M < 1.2)
///
/// Value: 0.015 (continuous with the supersonic base at M=1.2)
/// Physical basis: Shock waves disrupt circulation patterns, reducing Magnus effect
/// Effect: Spin drift significantly reduced in transonic flight
/// Source: Experimental spinning projectile studies
const MAGNUS_COEFF_TRANSONIC_REDUCTION: f64 = 0.015;

/// Base Magnus coefficient for supersonic flow (M > 1.2)
///
/// Value: 0.015 (dimensionless coefficient)
/// Physical basis: Shock-dominated flow with reduced but persistent circulation
/// Effect: Lower Magnus effect than subsonic, but higher than transonic minimum
const MAGNUS_COEFF_SUPERSONIC_BASE: f64 = 0.015;

/// Magnus coefficient scaling factor for high supersonic speeds
///
/// Value: 0.0044 (additional scaling with Mach number)
/// Formula: Magnus_coeff = BASE + SCALE * (M - 1.2) for M > 1.2
/// Physical basis: Partial recovery of circulation effects at higher Mach numbers
const MAGNUS_COEFF_SUPERSONIC_SCALE: f64 = 0.0044;

/// Transonic regime boundaries for Magnus effect calculations
const MAGNUS_TRANSONIC_LOWER: f64 = 0.8; // Lower bound of transonic regime
const MAGNUS_TRANSONIC_UPPER: f64 = 1.2; // Upper bound of transonic regime
const MAGNUS_TRANSONIC_RANGE: f64 = 0.4; // Range width (1.2 - 0.8)
const MAGNUS_SUPERSONIC_RANGE: f64 = 1.8; // Scaling range for supersonic recovery

// Note: These Magnus coefficients are calibrated against real-world spin drift measurements
// from McCoy's "Modern Exterior Ballistics" and experimental data. The dimensionless
// coefficients represent the Magnus moment per unit angle of attack.

// Atmosphere detection thresholds
const MAX_REALISTIC_DENSITY: f64 = 2.0; // kg/m³
const MIN_REALISTIC_SPEED_OF_SOUND: f64 = 200.0; // m/s

fn dry_air_temperature_c_from_sound_speed(speed_of_sound_mps: f64) -> f64 {
    speed_of_sound_mps * speed_of_sound_mps / (1.4 * 287.05) - 273.15
}

/// Calculate Magnus moment coefficient C_Lα based on Mach number
/// Based on McCoy's 'Modern Exterior Ballistics' and empirical data
pub(crate) fn calculate_magnus_moment_coefficient(mach: f64) -> f64 {
    // Magnus moment coefficient varies with Mach number
    // Values based on empirical data for spitzer bullets

    if mach < MAGNUS_TRANSONIC_LOWER {
        // Subsonic: relatively constant
        MAGNUS_COEFF_SUBSONIC
    } else if mach < MAGNUS_TRANSONIC_UPPER {
        // Transonic: reduced due to shock formation
        // Linear interpolation through transonic region
        MAGNUS_COEFF_SUBSONIC
            - MAGNUS_COEFF_TRANSONIC_REDUCTION * (mach - MAGNUS_TRANSONIC_LOWER)
                / MAGNUS_TRANSONIC_RANGE
    } else {
        // Supersonic: gradually recovers
        MAGNUS_COEFF_SUPERSONIC_BASE
            + MAGNUS_COEFF_SUPERSONIC_SCALE
                * ((mach - MAGNUS_TRANSONIC_UPPER) / MAGNUS_SUPERSONIC_RANGE).min(1.0)
    }
}

/// Project a vector expressed in level downrange/up/lateral axes into the inclined shot frame.
///
/// Shot-frame X follows the slant trajectory, Y is perpendicular to it in the vertical plane,
/// and Z remains lateral. The zero-angle return preserves level-fire values bit-for-bit.
#[inline]
pub(crate) fn level_vector_to_shot_frame(
    vector: Vector3<f64>,
    shooting_angle: f64,
) -> Vector3<f64> {
    if shooting_angle == 0.0 {
        return vector;
    }

    let (sin_angle, cos_angle) = shooting_angle.sin_cos();
    Vector3::new(
        vector.x * cos_angle + vector.y * sin_angle,
        -vector.x * sin_angle + vector.y * cos_angle,
        vector.z,
    )
}

/// Direction of the Magnus force generated by the lateral yaw of repose.
///
/// In the McCoy frame the yaw of repose is lateral, so `F_M ∝ v̂ × α_R` reduces to
/// gravity projected onto the plane normal to flight: down for right-hand twist and up for
/// left-hand twist. Passing the actual shot-frame gravity vector keeps this correct for inclined
/// fire; a vertical shot has no normal gravity component and therefore no repose-driven Magnus
/// direction.
pub(crate) fn yaw_of_repose_magnus_direction(
    air_velocity: Vector3<f64>,
    gravity_acceleration: Vector3<f64>,
    is_twist_right: bool,
) -> Option<Vector3<f64>> {
    let speed = air_velocity.norm();
    if !speed.is_finite() || speed <= 1e-12 {
        return None;
    }

    let velocity_unit = air_velocity / speed;
    let gravity_normal =
        gravity_acceleration - velocity_unit * gravity_acceleration.dot(&velocity_unit);
    let normal_magnitude = gravity_normal.norm();
    if !normal_magnitude.is_finite() || normal_magnitude <= 1e-12 {
        return None;
    }

    let right_twist_direction = gravity_normal / normal_magnitude;
    Some(if is_twist_right {
        right_twist_direction
    } else {
        -right_twist_direction
    })
}

/// Compute ballistic derivatives for trajectory integration.
///
/// `wind_vector` and `omega_vector` use level downrange/up/lateral axes and are projected into
/// the inclined shot frame internally.
#[allow(clippy::too_many_arguments)]
pub fn compute_derivatives(
    pos: Vector3<f64>,
    vel: Vector3<f64>,
    inputs: &BallisticInputs,
    wind_vector: Vector3<f64>,
    atmos_params: (f64, f64, f64, f64),
    bc_used: f64,
    omega_vector: Option<Vector3<f64>>,
    // MBA-1134: the in-integration spin-drift term that consumed `time` is deprecated (spin drift
    // is now a Litz post-process), so this is currently unused; kept in the signature for callers.
    _time: f64,
    // MBA-1137: optional downrange-segmented atmosphere. When `Some`, the STANDARD-mode base
    // (station-referenced) T/P/H is swapped for the zone selected by downrange distance (pos.x)
    // before the altitude lapse; the direct-atmosphere sentinel path is left untouched. `None`
    // (every existing caller) is byte-identical to the pre-feature behavior.
    atmo_sock: Option<&AtmoSock>,
) -> [f64; 6] {
    // Gravity acceleration vector, rotated into the shot-aligned frame by shooting_angle
    // (uphill/downhill inclined fire), matching cli_api::TrajectorySolver::gravity_acceleration.
    let theta = inputs.shooting_angle;
    let accel_gravity = Vector3::new(
        -G_ACCEL_MPS2 * theta.sin(),
        -G_ACCEL_MPS2 * theta.cos(),
        0.0,
    );

    // Wind-adjusted velocity. Wind sources are defined in the level frame; rotate the selected
    // vector into the same inclined shot frame used by velocity and gravity.
    let wind_vector = level_vector_to_shot_frame(wind_vector, theta);
    let velocity_adjusted = vel - wind_vector;
    let speed_air = velocity_adjusted.norm();

    // Initialize drag acceleration
    let mut accel_drag = Vector3::zeros();
    let mut accel_magnus = Vector3::zeros();

    // Calculate drag if velocity is significant
    if speed_air > crate::constants::MIN_VELOCITY_THRESHOLD {
        let v_rel_fps = speed_air * MPS_TO_FPS;

        // Get atmospheric conditions
        let altitude_at_pos = crate::atmosphere::shot_frame_altitude(
            inputs.altitude,
            pos[0],
            pos[1],
            inputs.shooting_angle,
        );

        // Check if we have direct atmosphere values
        // Direct atmosphere is indicated by having only 2 parameters where:
        // params[0] = air density, params[1] = speed of sound
        // params[2] and params[3] would be 0.0
        // BUT: we need to check if params[0] is a reasonable density value (< 2.0 kg/m³)
        let (air_density, speed_of_sound, temperature_c) = if atmos_params.0 < MAX_REALISTIC_DENSITY
            && atmos_params.1 > MIN_REALISTIC_SPEED_OF_SOUND
            && atmos_params.2 == 0.0
            && atmos_params.3 == 0.0
        {
            // Direct atmosphere values: atmos_params.1 is the SPEED OF SOUND here, NOT Celsius,
            // so back-compute temperature from it (c = sqrt(1.4*287.05*T_k)) for the optional
            // Reynolds-correction input below.
            let (rho, sound) = get_direct_atmosphere(atmos_params.0, atmos_params.1);
            (rho, sound, dry_air_temperature_c_from_sound_speed(sound))
        } else {
            // Calculate from base parameters. MBA-1136 (rank 9): route the local speed of sound
            // through the moist-air formula. Humidity is NOT plumbed to this call site — on the
            // sole production path (trajectory_integration::build_inputs) the `humidity` FIELD is
            // overwritten with atmos_params.3 (the density RATIO), so `inputs.humidity` here is not
            // a real RH. Pass 0.0 (dry) rather than fabricate humidity; density is unchanged and
            // the dry speed of sound is numerically identical to the old get_local_atmosphere.
            //
            // MBA-1137: when a downrange-segmented atmosphere is present, swap the BASE
            // (station-referenced) temp/pressure/ratio for the zone selected by downrange distance
            // (pos[0]) BEFORE the altitude lapse. The zone base_ratio is recomputed via CIPM from
            // the zone (temp, pressure, humidity); the swapped base then flows through the same
            // altitude-lapse pipeline, so downrange-zone selection and the world-vertical lapse
            // compose without double-counting.
            // Only the base density changes — the local speed of sound is independent of base_ratio,
            // so it still tracks the lapsed temperature/pressure. `None` -> byte-identical.
            let (base_temp_c, base_press_hpa, base_ratio) = match atmo_sock {
                Some(sock) => {
                    let (zt, zp, zh) = sock.atmo_for_range(pos[0]);
                    (zt, zp, calculate_air_density_cimp(zt, zp, zh) / 1.225)
                }
                None => {
                    // A zero/nonpositive standard-mode ratio means "not supplied", not vacuum.
                    // Match the plain fast path's sea-level fallback so sibling solvers cannot
                    // interpret the same FastIntegrationParams tuple oppositely (MBA-1157).
                    let base_ratio = if atmos_params.3 > 0.0 {
                        atmos_params.3
                    } else {
                        1.0
                    };
                    (atmos_params.1, atmos_params.2, base_ratio)
                }
            };
            let (rho, sound) = get_local_atmosphere_humid(
                altitude_at_pos,
                atmos_params.0, // base_alt
                base_temp_c,    // base_temp_c (zone-swapped when AtmoSock present)
                base_press_hpa, // base_press_hpa (zone-swapped when AtmoSock present)
                base_ratio,     // base_ratio (zone-swapped when AtmoSock present)
                0.0,            // humidity not available here (see note above)
            );
            // LOCAL temperature at the projectile altitude, back-computed from the LOCAL speed of
            // sound (get_local_atmosphere returns density/sound at altitude_at_pos but not temp;
            // its sound = sqrt(1.4*287.05*T_k)). Using base_temp_c here would feed the optional
            // Reynolds viscosity input a shooter-altitude temperature while density/sound are
            // local if that correction is enabled later.
            (rho, sound, dry_air_temperature_c_from_sound_speed(sound))
        };

        // Calculate Mach number with safe division
        let mach = if speed_of_sound > 1e-9 {
            speed_air / speed_of_sound
        } else {
            0.0 // No meaningful Mach number at zero speed of sound
        };

        // Get the base drag coefficient; optional corrections are applied or controlled below.
        let drag_factor = get_drag_coefficient_full(
            mach,
            &inputs.bc_type,
            false, // transonic applied exactly once below (was double-applied here + in block)
            false, // Reynolds remains disabled consistently across all solver families (MBA-945)
            None,  // let it determine shape
            if inputs.caliber_inches > 0.0 {
                Some(inputs.caliber_inches)
            } else {
                Some(inputs.bullet_diameter / 0.0254) // meters -> inches
            },
            if inputs.weight_grains > 0.0 {
                Some(inputs.weight_grains)
            } else {
                Some(inputs.bullet_mass / crate::constants::GRAINS_TO_KG) // kg -> grains
            },
            Some(speed_air),
            Some(air_density),
            Some(temperature_c),
        );

        // Get BC value
        let mut bc_val = bc_used;

        if inputs.use_bc_segments {
            // First try velocity-based segments if available
            if inputs.bc_segments_data.is_some() {
                bc_val = get_bc_for_velocity(v_rel_fps, inputs, bc_used);
            } else {
                match inputs.bc_segments.as_deref() {
                    // Fall back to a non-empty Mach table when no velocity data exist.
                    Some(segments) if !segments.is_empty() => {
                        bc_val = interpolated_bc(mach, segments, Some(inputs));
                    }
                    // An explicitly empty table is a no-op and preserves the active caller BC.
                    Some(_) => {}
                    // No explicit table: retain the opt-in automatic estimation behavior.
                    None => bc_val = get_bc_for_velocity(v_rel_fps, inputs, bc_used),
                }
            }
        } else if let Some(segments) = inputs
            .bc_segments
            .as_deref()
            .filter(|segments| !segments.is_empty())
        {
            // Explicit Mach-based segments (legacy behavior when use_bc_segments=false)
            bc_val = interpolated_bc(mach, segments, Some(inputs));
        }

        // Guard bc_val == 0 (allowed on the FFI/WASM/library surfaces, which lack the CLI's
        // 0.001 floor, and a user-supplied BC segment can be 0): the drag division below would be
        // Inf -> NaN, poisoning the whole trajectory. Mirrors the guards already in
        // cli_api::calculate_acceleration and fast_trajectory::compute_derivatives. Inert for
        // valid BCs (>= 0.001).
        let bc_val = bc_val.max(1e-6);

        // Apply the documented radian tip-off yaw, decaying exponentially with distance.
        let yaw_rad = if inputs.tipoff_decay_distance.abs() > 1e-9 {
            inputs.tipoff_yaw * (-pos[0] / inputs.tipoff_decay_distance).exp()
        } else {
            inputs.tipoff_yaw // No decay if distance is zero
        };
        // MBA-1227: quadratic yaw drag is ADDITIVE per McCoy: CD = CD0 + CD_delta2*delta^2.
        // The previous multiplicative form Cd*(1 + delta^2) implied CD_delta2 == CD0
        // (~0.3 per rad^2), an order of magnitude below literature (~4-20 per rad^2 for
        // spitzer rifle bullets). The additive term is expressed in retardation space
        // (actual-Cd over sectional density), which is denominator-correct for BOTH the
        // G-model path (drag_factor/BC == Cd_actual/SD) and the custom-table path
        // (drag_factor/SD). Computed after the (drag_factor, retard_denom) resolution below.

        // Calculate density scaling
        let density_scale = air_density / STANDARD_AIR_DENSITY;

        // Apply the transonic drag-rise correction exactly ONCE. The base Cd above is taken
        // WITHOUT transonic correction (apply_transonic_correction=false), so this is the only
        // application. Previously the correction was applied here AND inside
        // get_drag_coefficient_full, which squared the drag-rise factor and double-counted wave
        // drag across the transonic band (Cd ~3x too high near Mach 1). transonic_correction
        // self-gates via the projectile's critical Mach (returns the input unchanged outside the
        // band), and include_wave_drag=false matches cli_api::calculate_drag_coefficient — the
        // G1/G7 tables already embed the transonic rise, so additive wave drag would double-count.
        // Use the same SI fallbacks as the get_drag_coefficient_full call above (and
        // fast_trajectory): an SI-only caller may leave caliber_inches/weight_grains at 0, so
        // derive them from the SI bullet_diameter/bullet_mass rather than feeding zeros into
        // get_projectile_shape (which would mis-classify the shape via weight/caliber).
        let caliber_in = if inputs.caliber_inches > 0.0 {
            inputs.caliber_inches
        } else {
            inputs.bullet_diameter / 0.0254 // meters -> inches
        };
        let weight_gr = if inputs.weight_grains > 0.0 {
            inputs.weight_grains
        } else {
            inputs.bullet_mass / crate::constants::GRAINS_TO_KG // kg -> grains
        };
        // MBA-949: shared resolver so named bullet_model shapes are honored here too (this path
        // previously used only the caliber/weight heuristic and ignored the name).
        let shape = crate::transonic_drag::resolve_projectile_shape(
            inputs.bullet_model.as_deref(),
            caliber_in,
            weight_gr,
            &inputs.bc_type.to_string(),
        );
        let drag_factor =
            crate::transonic_drag::transonic_correction(mach, drag_factor, shape, false);

        // MBA-945: the low-velocity Reynolds drag correction was applied ONLY here (not in cli_api
        // or fast_trajectory), so subsonic shots diverged across the three solver families. Removed
        // for consistency — py_ballisticcalc (the validation reference) does not model it, and
        // cli_api already matches pbc subsonically without it (validated to 1300yd in MBA-939). The
        // reynolds module and get_drag_coefficient_full's apply_reynolds flag remain available for a
        // future opt-in wired across all three solvers.

        // MBA-940: a user-supplied custom drag table overrides the G-model Cd entirely and is used
        // as-is — no reference-table transonic correction or name-derived form-factor multiplier
        // is applied to it (the curve already encodes the projectile's true drag, so applying one
        // would distort/double-count it).
        // The custom table's Cd is the projectile's ACTUAL drag coefficient, so the
        // retardation denominator must be the sectional density (lb/in²), not a BC:
        // Cd_own / SD == Cd_ref / BC (see BallisticInputs::custom_drag_denominator).
        let (drag_factor, retard_denom) = match inputs.custom_drag_table {
            Some(ref table) => (
                // MBA-1357: cd_scale is a single whole-curve drag multiplier applied here, at
                // the Cd lookup site. The Mach-keyed DSF table (truing_dsf.rs) is a SEPARATE,
                // drop-only post-processing correction applied to a solved TrajectoryResult's
                // points after integration finishes — it never touches this drag computation.
                table.interpolate(mach) * inputs.cd_scale,
                inputs.custom_drag_denominator(bc_val),
            ),
            None => (drag_factor, bc_val),
        };

        // Calculate drag acceleration
        let standard_factor = drag_factor * CD_TO_RETARD;
        let mut a_drag_ft_s2 =
            (v_rel_fps.powi(2) * standard_factor * density_scale) / retard_denom;
        // MBA-1227 additive yaw-drag term (see comment above). Skipped entirely when
        // tip-off yaw is zero (the default), leaving baseline trajectories bit-identical.
        // When sectional density is unavailable (degenerate SI-less inputs) the term is
        // dropped rather than mis-scaled through an unrelated denominator.
        if yaw_rad != 0.0 {
            if let Some(sd) = inputs.sectional_density_lb_in2() {
                if sd > 0.0 {
                    a_drag_ft_s2 += v_rel_fps.powi(2)
                        * CD_TO_RETARD
                        * density_scale
                        * (inputs.cd_delta2 * yaw_rad.powi(2) / sd);
                }
            }
        }
        let a_drag_m_s2 = a_drag_ft_s2 * FPS_TO_MPS;

        // Apply drag in opposite direction of relative velocity
        accel_drag = -a_drag_m_s2 * (velocity_adjusted / speed_air);

        // Magnus Effect calculation. Gated on enable_magnus specifically so it is
        // independent of Coriolis (matches the cli_api solver's decoupled flags).
        // MBA-1134 (rank 35): when the canonical empirical Litz spin-drift post-process is active
        // (use_enhanced_spin_drift + advanced effects — the same condition that drove the now-
        // deprecated in-integration spin-drift term), it already captures the gyroscopic/yaw-of-
        // repose lateral. The explicit Magnus side force must NOT be added on top or the two
        // lateral models stack and double-count the drift, so suppress Magnus in that case.
        if inputs.enable_magnus
            && !(inputs.use_enhanced_spin_drift && inputs.enable_advanced_effects)
            && inputs.bullet_diameter > 0.0
            && inputs.twist_rate > 0.0
        {
            let diameter_m = inputs.bullet_diameter;
            let (spin_rate_rad_s, spin_param) = crate::spin_drift::calculate_magnus_spin_state(
                inputs.muzzle_velocity,
                speed_air,
                inputs.twist_rate,
                diameter_m,
            );

            let c_np = calculate_magnus_moment_coefficient(mach);

            // Calculate reference area
            let area = std::f64::consts::PI * (diameter_m / 2.0).powi(2);

            // Yaw of repose for the proper Magnus force. Stability/yaw helpers are
            // imperial: use the explicit imperial mirror fields, and convert the SI
            // bullet_length to inches at this boundary.
            let d_in = inputs.caliber_inches;
            let m_gr = inputs.weight_grains;
            let l_in = if inputs.bullet_length > 0.0 {
                inputs.bullet_length / 0.0254 // meters -> inches
            } else {
                // MBA-1135: mass-based length estimate (was a mass-blind 4.5-caliber default).
                let est_m = crate::stability::estimate_bullet_length_m(
                    inputs.bullet_diameter,
                    inputs.bullet_mass,
                );
                if est_m > 0.0 {
                    est_m / 0.0254
                } else {
                    4.5 * d_in.max(1e-9)
                }
            };
            // Use current-flight Sg with the muzzle-set spin. Back-calculating the effective
            // twist from fixed spin and current airspeed lets gyroscopic stability (and therefore
            // yaw of repose) grow as translational velocity decays; local density supplies the
            // canonical Miller atmospheric correction.
            let sg = crate::spin_drift::calculate_dynamic_stability(
                m_gr,
                speed_air,
                spin_rate_rad_s,
                d_in,
                l_in,
                air_density,
            );
            let (yaw_rad, _) = crate::spin_drift::calculate_yaw_of_repose(
                sg,
                speed_air,
                spin_rate_rad_s,
                0.0,
                0.0,
                air_density,
                d_in,
                l_in,
                m_gr,
                mach,
                "match",
                false,
            );

            // Proper McCoy Magnus FORCE: F = q S C_Npa (pd/2V) sin(alpha_R).
            let magnus_force_magnitude =
                0.5 * air_density * speed_air.powi(2) * area * c_np * spin_param * yaw_rad.sin();

            // A lateral yaw of repose produces a vertical Magnus force. Its lateral force is
            // aerodynamic lift, represented separately by the Litz spin-drift model.
            if magnus_force_magnitude > 1e-12 {
                if let Some(magnus_direction) = yaw_of_repose_magnus_direction(
                    velocity_adjusted,
                    accel_gravity,
                    inputs.is_twist_right,
                ) {
                    let bullet_mass_kg = inputs.bullet_mass; // already kg (SI)
                    accel_magnus = (magnus_force_magnitude / bullet_mass_kg) * magnus_direction;
                }
            }
        }
    }

    // Total acceleration
    let mut accel = accel_gravity + accel_drag + accel_magnus;

    // Add Coriolis acceleration if a level-frame omega vector is provided. The physical term is
    // -2 Ω×v (MBA-957: the old +2 "frame-relabel" justification was wrong — it flipped the
    // lateral drift; the caller now builds omega with the corrected lateral sign, matching the
    // validated cli_api solver, so the canonical -2 applies directly).
    if let Some(omega) = omega_vector {
        let omega = level_vector_to_shot_frame(omega, theta);
        let accel_coriolis = -2.0 * omega.cross(&vel);
        accel += accel_coriolis;
    }

    // MBA-1134 (rank 10): the in-integration enhanced-spin-drift ACCELERATION term
    // (spin_drift::calculate_enhanced_spin_drift / apply_enhanced_spin_drift) is DEPRECATED and is
    // no longer applied here. It carried a unit bug and was a SECOND, inconsistent lateral model.
    // Spin drift is now the single canonical empirical Litz post-process
    // (spin_drift::litz_drift_meters), applied by cli_api::apply_spin_drift and the fast /
    // Monte-Carlo path (fast_trajectory + monte_carlo) at the endpoint time-of-flight. Keeping only
    // one model here guarantees the three solver families agree on lateral drift and prevents the
    // Magnus + spin-drift double-count (see the Magnus gate above). The `calculate_enhanced_spin_drift`
    // / `apply_enhanced_spin_drift` helpers remain in spin_drift.rs for backward compatibility but
    // are no longer wired into any integration path.

    // Return state derivatives: [velocity, acceleration]
    [vel[0], vel[1], vel[2], accel[0], accel[1], accel[2]]
}

/// Interpolate a ballistic coefficient from Mach-keyed segments.
///
/// An empty optional table is not a request to invent a new BC: when `inputs` are present, the
/// caller's scalar `bc_value` is preserved. Standalone calls without inputs retain the historical
/// conservative fallback because no caller BC is available.
pub fn interpolated_bc(
    mach: f64,
    segments: &[(f64, f64)],
    inputs: Option<&BallisticInputs>,
) -> f64 {
    if segments.is_empty() {
        if let Some(inputs) = inputs {
            return inputs.bc_value;
        }
        return crate::constants::BC_FALLBACK_CONSERVATIVE;
    }

    if segments.len() == 1 {
        return segments[0].1;
    }

    // Ensure ascending-Mach order for interpolation. Fast path: when the segments are
    // already sorted (the common case — they are normalized once at construction), borrow
    // them and skip the per-call heap alloc + O(n log n) sort on the integration hot path.
    let sorted_segments: std::borrow::Cow<[(f64, f64)]> =
        if segments.windows(2).all(|w| w[0].0 <= w[1].0) {
            std::borrow::Cow::Borrowed(segments)
        } else {
            let mut v = segments.to_vec();
            v.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
            std::borrow::Cow::Owned(v)
        };

    // Handle out-of-range cases first
    if mach <= sorted_segments[0].0 {
        return sorted_segments[0].1;
    }
    if mach >= sorted_segments[sorted_segments.len() - 1].0 {
        return sorted_segments[sorted_segments.len() - 1].1;
    }

    // Find the appropriate segment using binary search
    let idx = sorted_segments.partition_point(|(m, _)| *m <= mach);
    if idx == 0 || idx >= sorted_segments.len() {
        // Should not happen given the checks above
        return sorted_segments[0].1;
    }

    let (mach1, bc1) = sorted_segments[idx - 1];
    let (mach2, bc2) = sorted_segments[idx];

    // Linear interpolation with safe division
    let denominator = mach2 - mach1;
    if denominator.abs() < crate::constants::MIN_DIVISION_THRESHOLD {
        return bc1; // Return first BC value if Mach values are identical
    }
    let t = (mach - mach1) / denominator;
    bc1 + t * (bc2 - bc1)
}

/// Get BC value for current velocity, supporting velocity-based BC segments
fn get_bc_for_velocity(velocity_fps: f64, inputs: &BallisticInputs, bc_used: f64) -> f64 {
    // Check if velocity-based BC segments are enabled
    if !inputs.use_bc_segments {
        return bc_used;
    }

    // Try direct BC segments data first
    if let Some(bc_segments_data) = inputs
        .bc_segments_data
        .as_ref()
        .filter(|segments| !segments.is_empty())
    {
        // An explicit table is authoritative even when this velocity lies in a coverage gap.
        // Do not silently replace it with an auto-estimated table after a miss.
        return velocity_segment_bc(velocity_fps, bc_segments_data, bc_used);
    }

    // Try BC estimation if we have bullet details but no segments. MBA-955: the estimation is
    // factored into estimate_bc_segments_for so the per-integration setup (build_inputs) can
    // pre-populate bc_segments_data ONCE rather than rebuilding it here every step.
    if let Some(segments) = estimate_bc_segments_for(inputs, bc_used) {
        return velocity_segment_bc(velocity_fps, &segments, bc_used);
    }

    // Fallback to constant BC
    bc_used
}

/// Estimate velocity-BC segments from bullet characteristics (MBA-955). Extracted from
/// get_bc_for_velocity's slow path so the per-integration setup can compute the segments ONCE
/// (build_inputs pre-populates bc_segments_data) instead of rebuilding them — allocating a model
/// String and a segment Vec — on every derivative evaluation. Returns None when the bullet
/// details needed for estimation are absent (the caller then falls back to the constant BC). The
/// logic is byte-identical to the previous inline slow path.
pub(crate) fn estimate_bc_segments_for(
    inputs: &BallisticInputs,
    bc_used: f64,
) -> Option<Vec<crate::BCSegmentData>> {
    if !(inputs.bullet_diameter > 0.0 && inputs.bullet_mass > 0.0 && bc_used > 0.0) {
        return None;
    }
    // Model string from bullet_id or a generic weight-based description (unchanged).
    let model = if let Some(ref bullet_id) = inputs.bullet_id {
        bullet_id.clone()
    } else {
        format!("{}gr bullet", inputs.weight_grains as i32)
    };
    // Prefer the legacy explicit string when present, but otherwise preserve the
    // typed drag model carried by every current BallisticInputs constructor.
    let bc_type_str = inputs.bc_type_str.as_deref().unwrap_or(match inputs.bc_type {
        crate::DragModel::G7 => "G7",
        _ => "G1",
    });
    Some(BCSegmentEstimator::estimate_bc_segments(
        bc_used,
        inputs.caliber_inches,
        inputs.weight_grains,
        &model,
        bc_type_str,
    ))
}

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

    fn expected_shot_frame_vector(level: Vector3<f64>, angle: f64) -> Vector3<f64> {
        let (sin_angle, cos_angle) = angle.sin_cos();
        Vector3::new(
            level.x * cos_angle + level.y * sin_angle,
            -level.x * sin_angle + level.y * cos_angle,
            level.z,
        )
    }

    #[test]
    fn level_vector_projection_is_bit_exact_at_zero_incline() {
        let level = Vector3::new(12.5, -0.0, -7.25);
        let projected = level_vector_to_shot_frame(level, 0.0);

        for component in 0..3 {
            assert_eq!(projected[component].to_bits(), level[component].to_bits());
        }
    }

    #[test]
    fn dry_air_sound_speed_round_trips_to_local_celsius() {
        for expected_temp_c in [-40.0_f64, 15.0, 40.0] {
            let sound_speed = (1.4 * 287.05 * (expected_temp_c + 273.15)).sqrt();
            let actual_temp_c = dry_air_temperature_c_from_sound_speed(sound_speed);
            assert!(
                (actual_temp_c - expected_temp_c).abs() < 1e-10,
                "sound speed {sound_speed} m/s recovered {actual_temp_c} C, expected {expected_temp_c} C"
            );
        }

        let direct_mode_temp_c = dry_air_temperature_c_from_sound_speed(340.0);
        assert!(
            direct_mode_temp_c > 10.0 && direct_mode_temp_c < 20.0,
            "direct-mode 340 m/s must be interpreted as sound speed, not 340 C"
        );
    }

    #[test]
    fn tipoff_yaw_uses_documented_radians_for_drag_and_decay() {
        let mut baseline_inputs = create_test_inputs();
        baseline_inputs.tipoff_yaw = 0.0;
        baseline_inputs.tipoff_decay_distance = 50.0;
        let mut yawed_inputs = baseline_inputs.clone();
        yawed_inputs.tipoff_yaw = 0.1; // radians, per the field docs

        let drag_x = |inputs: &BallisticInputs, downrange_m: f64| {
            compute_derivatives(
                Vector3::new(downrange_m, 0.0, 0.0),
                Vector3::new(300.0, 0.0, 0.0),
                inputs,
                Vector3::zeros(),
                (1.225, 340.0, 0.0, 0.0),
                0.5,
                None,
                0.0,
                None,
            )[3]
        };

        // MBA-1227: yaw drag is ADDITIVE per McCoy (CD = CD0 + CD_delta2 * delta^2), so the
        // extra retardation has a closed form we can assert exactly:
        //   extra_mps2 = v_fps^2 * CD_TO_RETARD * density_scale * CD_delta2 * delta_eff^2 / SD
        // with delta_eff = tipoff_yaw * exp(-x / decay). At x = decay*ln(2) the effective yaw
        // halves, so the extra drag falls to a quarter — which also proves the documented
        // RADIAN interpretation and the decay wiring in one pass.
        let sd = baseline_inputs
            .sectional_density_lb_in2()
            .expect("test inputs carry mass/diameter");
        let v_fps = 300.0 * MPS_TO_FPS; // same (rounded) constant the kernel uses
        let density_scale = 1.225 / STANDARD_AIR_DENSITY;

        for (downrange_m, decay_scale) in [(0.0_f64, 1.0_f64), (50.0 * std::f64::consts::LN_2, 0.5)]
        {
            let delta_eff = 0.1 * decay_scale;
            let expected_extra_mps2 = v_fps.powi(2)
                * CD_TO_RETARD
                * density_scale
                * (yawed_inputs.cd_delta2 * delta_eff.powi(2) / sd)
                * FPS_TO_MPS;

            let baseline_drag = drag_x(&baseline_inputs, downrange_m);
            let yawed_drag = drag_x(&yawed_inputs, downrange_m);
            assert!(baseline_drag < 0.0, "baseline must be downrange drag");
            let measured_extra = baseline_drag - yawed_drag; // yawed is MORE negative

            assert!(
                (measured_extra - expected_extra_mps2).abs() <= 1e-9 * expected_extra_mps2,
                "additive yaw drag at x={downrange_m} m: measured extra {measured_extra:.9} \
                 != expected {expected_extra_mps2:.9}"
            );
        }
    }

    fn create_test_inputs() -> BallisticInputs {
        // SI-canonical geometry/mass (kg, meters) — same convention as the struct
        // docs and cli_api — plus the explicit imperial mirror fields
        // (caliber_inches/weight_grains) the stability/Magnus helpers read.
        BallisticInputs {
            muzzle_velocity: 800.0, // m/s
            bc_value: 0.5,
            bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG, // kg (168 gr)
            bullet_diameter: 0.308 * 0.0254,    // meters (.308 in)
            bullet_length: 1.215 * 0.0254,      // meters
            caliber_inches: 0.308,
            weight_grains: 168.0,
            altitude: 1000.0,
            ..Default::default()
        }
    }

    #[test]
    fn measured_bc_drag_ignores_name_based_form_factor_flag() {
        let derivatives_with_flag = |use_form_factor| {
            let inputs = BallisticInputs {
                bc_value: 0.462,
                bc_type: crate::DragModel::G1,
                bullet_model: Some("168gr SMK Match".to_string()),
                use_form_factor,
                ..create_test_inputs()
            };

            compute_derivatives(
                Vector3::zeros(),
                Vector3::new(600.0, 0.0, 0.0),
                &inputs,
                Vector3::zeros(),
                (1.225, 340.0, 0.0, 0.0),
                inputs.bc_value,
                None,
                0.0,
                None,
            )
        };

        let baseline = derivatives_with_flag(false);
        let flagged = derivatives_with_flag(true);

        for component in 3..6 {
            assert_eq!(
                flagged[component].to_bits(),
                baseline[component].to_bits(),
                "published BC already encodes form factor: component {component}, baseline={} flagged={}",
                baseline[component],
                flagged[component]
            );
        }
    }

    #[test]
    fn inclined_positions_at_same_world_altitude_have_same_atmospheric_acceleration() {
        let angle = std::f64::consts::FRAC_PI_6;
        let mut inputs = create_test_inputs();
        inputs.altitude = 100.0;
        inputs.shooting_angle = angle;
        let velocity = Vector3::new(600.0, 0.0, 0.0);
        let along_slant = Vector3::new(1_000.0, 0.0, 0.0);
        let across_slant = Vector3::new(0.0, 500.0 / angle.cos(), 0.0);
        let atmo = (inputs.altitude, 15.0, 1013.25, 1.0);

        let a = compute_derivatives(
            along_slant,
            velocity,
            &inputs,
            Vector3::zeros(),
            atmo,
            inputs.bc_value,
            None,
            0.0,
            None,
        );
        let b = compute_derivatives(
            across_slant,
            velocity,
            &inputs,
            Vector3::zeros(),
            atmo,
            inputs.bc_value,
            None,
            0.0,
            None,
        );

        for component in 3..6 {
            assert!(
                (a[component] - b[component]).abs() < 1e-10,
                "derivative component {component} differs at equal world altitude: {} vs {}",
                a[component],
                b[component]
            );
        }
    }

    #[test]
    fn inclined_headwind_is_rotated_into_solver_frame() {
        let angle = std::f64::consts::FRAC_PI_6;
        let mut inputs = create_test_inputs();
        inputs.shooting_angle = angle;
        let level_headwind = Vector3::new(-100.0, 0.0, 0.0);
        let velocity = expected_shot_frame_vector(level_headwind, angle);
        let actual = compute_derivatives(
            Vector3::zeros(),
            velocity,
            &inputs,
            level_headwind,
            (1.225, 340.0, 0.0, 0.0),
            inputs.bc_value,
            None,
            0.0,
            None,
        );
        let expected = Vector3::new(
            -G_ACCEL_MPS2 * angle.sin(),
            -G_ACCEL_MPS2 * angle.cos(),
            0.0,
        );

        assert!(
            (Vector3::new(actual[3], actual[4], actual[5]) - expected).norm() < 1e-12,
            "co-moving horizontal wind must leave only shot-frame gravity: {actual:?}"
        );
    }

    #[test]
    fn inclined_coriolis_is_rotated_into_solver_frame() {
        let angle = std::f64::consts::FRAC_PI_6;
        let mut inputs = create_test_inputs();
        inputs.shooting_angle = angle;
        let velocity = Vector3::new(600.0, 20.0, 5.0);
        let level_omega = Vector3::new(3.0e-5, 6.0e-5, -2.0e-5);
        let run = |omega| {
            compute_derivatives(
                Vector3::zeros(),
                velocity,
                &inputs,
                Vector3::zeros(),
                (1.225, 340.0, 0.0, 0.0),
                inputs.bc_value,
                omega,
                0.0,
                None,
            )
        };
        let baseline = run(None);
        let with_coriolis = run(Some(level_omega));
        let actual = Vector3::new(
            with_coriolis[3] - baseline[3],
            with_coriolis[4] - baseline[4],
            with_coriolis[5] - baseline[5],
        );
        let expected = -2.0 * expected_shot_frame_vector(level_omega, angle).cross(&velocity);

        assert!(
            (actual - expected).norm() < 1e-12,
            "inclined Coriolis mismatch: actual={actual:?}, expected={expected:?}"
        );
    }

    #[test]
    fn explicit_velocity_segment_gap_uses_scalar_bc_without_auto_estimation() {
        let mut inputs = create_test_inputs();
        inputs.bc_value = 0.91;
        inputs.use_bc_segments = true;
        inputs.bc_segments_data = Some(vec![
            crate::BCSegmentData {
                velocity_min: 0.0,
                velocity_max: 999.0,
                bc_value: 0.2,
            },
            crate::BCSegmentData {
                velocity_min: 1000.0,
                velocity_max: 2000.0,
                bc_value: 0.3,
            },
        ]);

        assert_eq!(
            get_bc_for_velocity(999.5, &inputs, inputs.bc_value),
            inputs.bc_value,
            "an explicit table gap must not be replaced by an auto-estimated segment"
        );
    }

    #[test]
    fn test_mba955_bc_segments_prepopulate_byte_identical() {
        // MBA-955: pre-populating bc_segments_data once (in build_inputs) must return
        // BYTE-IDENTICAL BC to the old per-step estimation. Build the slow-path inputs
        // (bc_segments_data = None -> get_bc_for_velocity estimates every call) and the
        // pre-populated inputs (bc_segments_data = estimate_bc_segments_for, the same helper
        // build_inputs now calls), and assert get_bc_for_velocity agrees bit-for-bit across the
        // whole velocity range.
        let mut slow = create_test_inputs();
        slow.use_bc_segments = true;
        slow.bc_segments_data = None;
        slow.bc_segments = None;

        let bc_used = slow.bc_value;
        let mut fast = slow.clone();
        fast.bc_segments_data = estimate_bc_segments_for(&fast, bc_used);
        assert!(
            fast.bc_segments_data.is_some(),
            "estimation should yield segments for a valid bullet"
        );

        for v in (200..=3500).step_by(50) {
            let vf = v as f64;
            let a = get_bc_for_velocity(vf, &slow, bc_used);
            let b = get_bc_for_velocity(vf, &fast, bc_used);
            assert_eq!(
                a.to_bits(),
                b.to_bits(),
                "BC differs at {vf} fps: slow={a} fast={b}"
            );
        }
    }

    #[test]
    fn estimated_segments_inherit_typed_g7_drag_model() {
        let mut inputs = create_test_inputs();
        inputs.bc_value = 0.243;
        inputs.bc_type = crate::DragModel::G7;
        inputs.bc_type_str = None;
        inputs.bullet_mass = 175.0 * crate::constants::GRAINS_TO_KG;
        inputs.weight_grains = 175.0;

        let actual = estimate_bc_segments_for(&inputs, inputs.bc_value).unwrap();
        let expected = BCSegmentEstimator::estimate_bc_segments(
            inputs.bc_value,
            inputs.caliber_inches,
            inputs.weight_grains,
            "175gr bullet",
            "G7",
        );

        assert_eq!(actual.len(), expected.len());
        for (actual, expected) in actual.iter().zip(&expected) {
            assert_eq!(actual.velocity_min.to_bits(), expected.velocity_min.to_bits());
            assert_eq!(actual.velocity_max.to_bits(), expected.velocity_max.to_bits());
            assert_eq!(actual.bc_value.to_bits(), expected.bc_value.to_bits());
        }

        // An explicitly populated legacy string remains authoritative.
        inputs.bc_type_str = Some("G1".to_string());
        let legacy = estimate_bc_segments_for(&inputs, inputs.bc_value).unwrap();
        let expected_g1 = BCSegmentEstimator::estimate_bc_segments(
            inputs.bc_value,
            inputs.caliber_inches,
            inputs.weight_grains,
            "175gr bullet",
            "G1",
        );
        assert_eq!(legacy.len(), expected_g1.len());
        for (legacy, expected) in legacy.iter().zip(&expected_g1) {
            assert_eq!(legacy.bc_value.to_bits(), expected.bc_value.to_bits());
        }
    }

    #[test]
    fn test_compute_derivatives_basic() {
        let pos = Vector3::new(0.0, 0.0, 0.0);
        let vel = Vector3::new(800.0, 0.0, 0.0);
        let inputs = create_test_inputs();
        let wind_vector = Vector3::zeros();
        // Use direct atmosphere values: (air_density, speed_of_sound, 0.0, 0.0)
        let atmos_params = (1.225, 340.0, 0.0, 0.0); // Standard air density and speed of sound
        let bc_used = 0.5;

        let result = compute_derivatives(
            pos,
            vel,
            &inputs,
            wind_vector,
            atmos_params,
            bc_used,
            None,
            0.0,
            None,
        );

        // Check that we get velocity and acceleration components
        assert_eq!(result.len(), 6);

        // Velocity components should match input velocity
        assert!((result[0] - vel[0]).abs() < 1e-10);
        assert!((result[1] - vel[1]).abs() < 1e-10);
        assert!((result[2] - vel[2]).abs() < 1e-10);

        // Should have gravitational acceleration
        assert!(result[4] < 0.0); // Negative y acceleration due to gravity

        // Should have drag acceleration opposing motion
        assert!(result[3] < 0.0); // Negative x acceleration due to drag
    }

    #[test]
    fn standard_atmosphere_zero_ratio_uses_sea_level_fallback() {
        let pos = Vector3::new(0.0, 0.0, 0.0);
        let vel = Vector3::new(800.0, 0.0, 0.0);
        let inputs = create_test_inputs();
        let run = |base_ratio| {
            compute_derivatives(
                pos,
                vel,
                &inputs,
                Vector3::zeros(),
                (inputs.altitude, 15.0, 1013.25, base_ratio),
                inputs.bc_value,
                None,
                0.0,
                None,
            )
        };

        let missing_ratio = run(0.0);
        let explicit_sea_level = run(1.0);
        assert!(
            missing_ratio[3] < 0.0,
            "missing standard density ratio must not produce vacuum drag"
        );
        assert!((missing_ratio[3] - explicit_sea_level[3]).abs() < 1e-12);
    }

    #[test]
    fn test_compute_derivatives_with_wind() {
        let pos = Vector3::new(0.0, 0.0, 0.0);
        let vel = Vector3::new(800.0, 0.0, 0.0);
        let inputs = create_test_inputs();
        let wind_vector = Vector3::new(10.0, 0.0, 0.0); // Tailwind
        let atmos_params = (1.225, 340.0, 0.0, 0.0); // Standard air density and speed of sound
        let bc_used = 0.5;

        let result = compute_derivatives(
            pos,
            vel,
            &inputs,
            wind_vector,
            atmos_params,
            bc_used,
            None,
            0.0,
            None,
        );

        // With tailwind, effective velocity should be lower, thus less drag
        // Just check that we have some drag (negative acceleration)
        assert!(result[3] < 0.0); // Should have drag
    }

    #[test]
    fn test_compute_derivatives_with_coriolis() {
        let pos = Vector3::new(0.0, 0.0, 0.0);
        let vel = Vector3::new(800.0, 0.0, 0.0);
        let inputs = create_test_inputs();
        let wind_vector = Vector3::zeros();
        let atmos_params = (1.225, 340.0, 0.0, 0.0); // Standard air density and speed of sound
        let bc_used = 0.5;
        let omega = Vector3::new(0.0, 0.0, 7.2921e-5); // Earth's rotation

        let result = compute_derivatives(
            pos,
            vel,
            &inputs,
            wind_vector,
            atmos_params,
            bc_used,
            Some(omega),
            0.0,
            None,
        );

        // Should have Coriolis effect
        assert!(result[4].abs() > 1e-3); // Should have some y-component from Coriolis
    }

    #[test]
    fn test_interpolated_bc() {
        let segments = vec![(0.5, 0.4), (1.0, 0.5), (1.5, 0.6), (2.0, 0.5)];

        // Test exact matches
        assert!((interpolated_bc(1.0, &segments, None) - 0.5).abs() < 1e-10);

        // Test interpolation
        let bc_075 = interpolated_bc(0.75, &segments, None);
        assert!(bc_075 > 0.4 && bc_075 < 0.5);

        // Test out of range
        assert!((interpolated_bc(0.1, &segments, None) - 0.4).abs() < 1e-10);
        assert!((interpolated_bc(3.0, &segments, None) - 0.5).abs() < 1e-10);
    }

    #[test]
    fn test_interpolated_bc_edge_cases() {
        // Empty segments
        assert!(
            (interpolated_bc(1.0, &[], None) - crate::constants::BC_FALLBACK_CONSERVATIVE).abs()
                < 1e-10
        );

        let mut inputs = create_test_inputs();
        inputs.bc_type = crate::DragModel::G7;
        inputs.bc_value = 0.487;
        assert_eq!(
            interpolated_bc(1.0, &[], Some(&inputs)).to_bits(),
            inputs.bc_value.to_bits(),
            "an empty optional table must preserve the caller's scalar BC"
        );

        // Single segment
        let single = vec![(1.0, 0.7)];
        assert!((interpolated_bc(1.5, &single, None) - 0.7).abs() < 1e-10);
    }

    #[test]
    fn empty_mach_segments_preserve_active_bc_used() {
        let mut inputs = create_test_inputs();
        inputs.bc_type = crate::DragModel::G7;
        inputs.bc_value = 0.123; // Deliberately different from the active fitted/adjusted BC.

        let drag_acceleration = |inputs: &BallisticInputs| {
            compute_derivatives(
                Vector3::zeros(),
                Vector3::new(700.0, 0.0, 0.0),
                inputs,
                Vector3::zeros(),
                (1.225, 340.0, 0.0, 0.0),
                0.487,
                None,
                0.0,
                None,
            )[3]
        };

        inputs.use_bc_segments = false;
        inputs.bc_segments = None;
        let no_table = drag_acceleration(&inputs);

        for use_bc_segments in [false, true] {
            inputs.use_bc_segments = use_bc_segments;
            inputs.bc_segments = Some(Vec::new());
            let empty_table = drag_acceleration(&inputs);

            assert_eq!(
                empty_table.to_bits(),
                no_table.to_bits(),
                "Some(empty) must preserve bc_used when use_bc_segments={use_bc_segments}"
            );
        }
    }

    #[test]
    fn test_magnus_effect() {
        let pos = Vector3::new(0.0, 0.0, 0.0);
        let vel = Vector3::new(822.96, 0.0, 0.0); // 2700 fps
        let wind_vector = Vector3::zeros();
        let atmos_params = (1.225, 340.0, 0.0, 0.0); // Standard air density and speed of sound
        let bc_used = 0.5;
        let acceleration = |enable_magnus, is_twist_right| {
            let mut inputs = create_test_inputs();
            inputs.twist_rate = 10.0; // 1:10 twist
            inputs.is_twist_right = is_twist_right;
            inputs.enable_magnus = enable_magnus; // decoupled from enable_advanced_effects

            let result = compute_derivatives(
                pos,
                vel,
                &inputs,
                wind_vector,
                atmos_params,
                bc_used,
                None,
                0.0,
                None,
            );
            Vector3::new(result[3], result[4], result[5])
        };

        let baseline = acceleration(false, true);
        let right_twist = acceleration(true, true) - baseline;
        let left_twist = acceleration(true, false) - baseline;

        // Yaw of repose is lateral, so its Magnus force is vertical: down for right-hand twist
        // and up for left-hand twist. The lateral force from yaw of repose is lift, represented by
        // the separate Litz spin-drift model when that model is enabled.
        assert!(
            right_twist.y < 0.0,
            "right-hand Magnus must point down, got {right_twist:?}"
        );
        assert!(
            left_twist.y > 0.0,
            "left-hand Magnus must point up, got {left_twist:?}"
        );
        assert!((right_twist.y + left_twist.y).abs() < 1e-12);
        assert!(right_twist.x.abs() < 1e-12 && right_twist.z.abs() < 1e-12);
        assert!(left_twist.x.abs() < 1e-12 && left_twist.z.abs() < 1e-12);
        assert!(
            right_twist.y.abs() < 0.05,
            "Magnus must remain a small force"
        );
    }

    #[test]
    fn magnus_uses_velocity_corrected_muzzle_stability_gate() {
        let muzzle_velocity = 1_400.0 / MPS_TO_FPS;
        let mut inputs = create_test_inputs();
        inputs.muzzle_velocity = muzzle_velocity;
        inputs.twist_rate = 15.0;
        inputs.enable_magnus = true;

        let bare_sg = crate::spin_drift::miller_stability(0.308, 168.0, 15.0, 1.215);
        let canonical_sg = crate::spin_drift::effective_sg_from_inputs(&inputs, 15.0, 1013.25);
        assert!(bare_sg > 1.0, "test requires bare Sg above the Magnus gate");
        assert!(
            canonical_sg < 1.0,
            "velocity-corrected Sg must be below the gate, got {canonical_sg}"
        );

        let acceleration = |inputs: &BallisticInputs| {
            let result = compute_derivatives(
                Vector3::zeros(),
                Vector3::new(muzzle_velocity, 0.0, 0.0),
                inputs,
                Vector3::zeros(),
                (1.225, 340.0, 0.0, 0.0),
                0.5,
                None,
                0.0,
                None,
            );
            Vector3::new(result[3], result[4], result[5])
        };
        let enabled = acceleration(&inputs);
        inputs.enable_magnus = false;
        let disabled = acceleration(&inputs);

        assert_eq!(
            enabled, disabled,
            "canonical Sg below 1 must suppress every Magnus acceleration component"
        );
    }

    #[test]
    fn magnus_force_grows_as_fixed_spin_projectile_slows() {
        let mut inputs = create_test_inputs();
        inputs.muzzle_velocity = 800.0;
        inputs.twist_rate = 12.0;
        inputs.enable_magnus = true;

        let magnus_acceleration = |speed_mps| {
            let evaluate = |enable_magnus| {
                let mut run_inputs = inputs.clone();
                run_inputs.enable_magnus = enable_magnus;
                compute_derivatives(
                    Vector3::zeros(),
                    Vector3::new(speed_mps, 0.0, 0.0),
                    &run_inputs,
                    Vector3::zeros(),
                    (1.225, 340.0, 0.0, 0.0),
                    0.5,
                    None,
                    0.0,
                    None,
                )[4]
            };
            (evaluate(true) - evaluate(false)).abs()
        };

        let fast = magnus_acceleration(200.0);
        let slow = magnus_acceleration(100.0);
        let ratio = slow / fast;
        let expected_ratio = 2.0_f64.powf(5.0 / 3.0);

        assert!(fast > 0.0 && slow > 0.0, "fast={fast}, slow={slow}");
        assert!(
            (ratio - expected_ratio).abs() < 1e-3,
            "fixed-spin Magnus acceleration must grow downrange; slow/fast={ratio}, \
             expected={expected_ratio}"
        );
    }

    #[test]
    fn test_magnus_moment_coefficient() {
        // Test at various Mach numbers with corrected coefficients
        assert!((calculate_magnus_moment_coefficient(0.5) - 0.030).abs() < 0.001); // Subsonic
        assert!((calculate_magnus_moment_coefficient(0.8) - 0.030).abs() < 0.001); // Start of transonic
        assert!((calculate_magnus_moment_coefficient(1.0) - 0.0225).abs() < 0.001); // Mid transonic
        assert!((calculate_magnus_moment_coefficient(1.2) - 0.015).abs() < 0.001); // End of transonic
        assert!((calculate_magnus_moment_coefficient(2.0) - 0.01653).abs() < 0.001);
        // Supersonic
    }

    /// MBA-1356: cd_scale — this module's custom-deck interpolation site.
    fn deck_test_inputs(cd_scale: f64) -> BallisticInputs {
        BallisticInputs {
            custom_drag_table: Some(crate::drag::DragTable::new(
                vec![0.5, 1.0, 2.0, 3.0],
                vec![0.23, 0.40, 0.30, 0.26],
            )),
            cd_scale,
            ..create_test_inputs()
        }
    }

    fn deck_accel_x(cd_scale: f64) -> f64 {
        let inputs = deck_test_inputs(cd_scale);
        compute_derivatives(
            Vector3::zeros(),
            Vector3::new(700.0, 0.0, 0.0),
            &inputs,
            Vector3::zeros(),
            (1.225, 340.0, 0.0, 0.0),
            inputs.bc_value,
            None,
            0.0,
            None,
        )[3]
    }

    #[test]
    fn cd_scale_default_is_one_and_absent_matches_explicit() {
        assert_eq!(BallisticInputs::default().cd_scale, 1.0);

        let omitted = BallisticInputs {
            custom_drag_table: Some(crate::drag::DragTable::new(
                vec![0.5, 1.0, 2.0, 3.0],
                vec![0.23, 0.40, 0.30, 0.26],
            )),
            ..create_test_inputs()
        };
        assert_eq!(omitted.cd_scale, 1.0);

        let a_omitted = compute_derivatives(
            Vector3::zeros(),
            Vector3::new(700.0, 0.0, 0.0),
            &omitted,
            Vector3::zeros(),
            (1.225, 340.0, 0.0, 0.0),
            omitted.bc_value,
            None,
            0.0,
            None,
        );
        let a_explicit = deck_accel_x(1.0);
        assert_eq!(
            a_omitted[3].to_bits(),
            a_explicit.to_bits(),
            "omitted cd_scale (Default) must be bit-identical to an explicit 1.0"
        );
    }

    /// (b) Scale-direction test, derivatives-driven path: cd_scale=1.10 must add more drag
    /// (a more negative x-acceleration) than 1.0; 0.90 must add less.
    #[test]
    fn cd_scale_direction_on_derivatives_kernel() {
        let baseline = deck_accel_x(1.0);
        let scaled_up = deck_accel_x(1.10);
        let scaled_down = deck_accel_x(0.90);

        assert!(
            scaled_up < baseline,
            "cd_scale=1.10 must increase drag deceleration (more negative ax): \
             base={baseline} up={scaled_up}"
        );
        assert!(
            scaled_down > baseline,
            "cd_scale=0.90 must decrease drag deceleration (less negative ax): \
             base={baseline} down={scaled_down}"
        );
    }

    /// cd_scale must be inert on the standard G-model/BC path (no custom_drag_table).
    #[test]
    fn cd_scale_is_inert_without_a_custom_drag_table() {
        let make = |cd_scale: f64| BallisticInputs {
            cd_scale,
            ..create_test_inputs()
        };
        let accel = |cd_scale: f64| {
            let inputs = make(cd_scale);
            compute_derivatives(
                Vector3::zeros(),
                Vector3::new(700.0, 0.0, 0.0),
                &inputs,
                Vector3::zeros(),
                (1.225, 340.0, 0.0, 0.0),
                inputs.bc_value,
                None,
                0.0,
                None,
            )[3]
        };
        assert_eq!(
            accel(1.0).to_bits(),
            accel(1.5).to_bits(),
            "cd_scale must not affect the G-model/BC drag path"
        );
    }
}