ballistics_engine/

spin_drift_advanced.rs

// Advanced spin drift model based on modern ballistics research
// Incorporates multiple empirical models from:
// - Bryan Litz's Applied Ballistics for Long Range Shooting
// - McCoy's Modern Exterior Ballistics
// - Courtney & Courtney spin drift research papers

use std::f64::consts::PI;

/// Advanced spin drift coefficients based on extensive field data
#[derive(Debug, Clone)]
pub struct SpinDriftCoefficients {
    /// Litz coefficient for gyroscopic drift (typically 0.8-1.5)
    pub litz_coefficient: f64,
    /// McCoy's aerodynamic jump factor
    pub mccoy_jump_factor: f64,
    /// Courtney's transonic adjustment
    pub transonic_factor: f64,
    /// Yaw damping coefficient
    pub yaw_damping: f64,
}

impl SpinDriftCoefficients {
    /// Get coefficients for specific bullet types based on empirical data
    pub fn for_bullet_type(bullet_type: &str) -> Self {
        match bullet_type.to_lowercase().as_str() {
            "match" | "bthp" | "boat_tail" => Self {
                litz_coefficient: 1.25,
                mccoy_jump_factor: 0.85,
                transonic_factor: 0.75,
                yaw_damping: 0.92,
            },
            "vld" | "very_low_drag" => Self {
                litz_coefficient: 1.15,
                mccoy_jump_factor: 0.78,
                transonic_factor: 0.68,
                yaw_damping: 0.88,
            },
            "hybrid" | "hybrid_ogive" => Self {
                litz_coefficient: 1.20,
                mccoy_jump_factor: 0.82,
                transonic_factor: 0.72,
                yaw_damping: 0.90,
            },
            "flat_base" | "fb" => Self {
                litz_coefficient: 1.35,
                mccoy_jump_factor: 0.95,
                transonic_factor: 0.85,
                yaw_damping: 0.95,
            },
            _ => Self::default(),
        }
    }

    pub fn default() -> Self {
        Self {
            litz_coefficient: 1.25,
            mccoy_jump_factor: 0.85,
            transonic_factor: 0.75,
            yaw_damping: 0.92,
        }
    }
}

/// Calculate advanced spin drift using multiple empirical models
pub fn calculate_advanced_spin_drift(
    stability_factor: f64,
    time_of_flight_s: f64,
    velocity_mps: f64,
    muzzle_velocity_mps: f64,
    spin_rate_rad_s: f64,
    caliber_m: f64,
    mass_kg: f64,
    air_density_kg_m3: f64,
    is_right_twist: bool,
    bullet_type: &str,
) -> f64 {
    // Edge cases: no drift if no time or no stability
    // MBA-205: Guard against division by zero
    if time_of_flight_s <= 0.0 || stability_factor <= 0.0
        || muzzle_velocity_mps <= 0.0 || air_density_kg_m3 <= 0.0 {
        return 0.0;
    }
    
    let coeffs = SpinDriftCoefficients::for_bullet_type(bullet_type);
    
    // Direction based on twist
    let sign = if is_right_twist { 1.0 } else { -1.0 };
    
    // Calculate Mach numbers
    let mach_current = velocity_mps / 343.0;
    let mach_muzzle = muzzle_velocity_mps / 343.0;
    
    // 1. Litz's empirical formula (primary component)
    let litz_drift = calculate_litz_drift(
        stability_factor,
        time_of_flight_s,
        coeffs.litz_coefficient,
    );
    
    // 2. McCoy's aerodynamic jump correction
    let jump_correction = calculate_aerodynamic_jump_correction(
        mach_muzzle,
        spin_rate_rad_s,
        caliber_m,
        mass_kg,
        coeffs.mccoy_jump_factor,
    );
    
    // 3. Transonic correction factor
    let transonic_correction = calculate_transonic_correction(
        mach_current,
        coeffs.transonic_factor,
    );
    
    // 4. Yaw damping effect
    let yaw_factor = calculate_yaw_damping_factor(
        stability_factor,
        time_of_flight_s,
        coeffs.yaw_damping,
    );
    
    // 5. Velocity decay correction (new research)
    let velocity_ratio = velocity_mps / muzzle_velocity_mps;
    let velocity_correction = velocity_ratio.powf(0.3);
    
    // Combine all components
    let total_drift = sign * (
        litz_drift * transonic_correction * yaw_factor * velocity_correction
        + jump_correction
    );
    
    // Apply atmospheric density correction
    let density_correction = (1.225 / air_density_kg_m3).sqrt();
    
    total_drift * density_correction
}

/// Litz's empirical formula with refined coefficients
fn calculate_litz_drift(stability: f64, time_s: f64, coefficient: f64) -> f64 {
    if stability <= 1.0 || time_s <= 0.0 {
        return 0.0;
    }
    
    // Refined Litz formula based on extensive field testing
    // SD = k * (SG + 1.2) * TOF^1.83
    // where k is empirically determined for bullet type
    let base_drift = coefficient * (stability + 1.2) * time_s.powf(1.83);
    
    // Convert inches to meters
    base_drift * 0.0254
}

/// McCoy's aerodynamic jump correction at muzzle exit
fn calculate_aerodynamic_jump_correction(
    mach: f64,
    spin_rate_rad_s: f64,
    caliber_m: f64,
    mass_kg: f64,
    jump_factor: f64,
) -> f64 {
    // Aerodynamic jump contributes to initial displacement
    // Based on McCoy's research on muzzle exit effects

    // MBA-205: Guard against division by zero when mach == 0
    if mach <= 0.0 {
        return 0.0;
    }

    let spin_parameter = spin_rate_rad_s * caliber_m / (2.0 * 343.0 * mach);
    
    // Jump magnitude in milliradians
    let jump_mrad = jump_factor * spin_parameter * (mass_kg / 0.01).sqrt();
    
    // Convert to lateral displacement (approximation for small angles)
    // This is a one-time displacement, not time-dependent
    jump_mrad * 0.001 * 100.0 // Approximate 100m reference distance
}

/// Transonic correction based on Courtney & Courtney research
fn calculate_transonic_correction(mach: f64, transonic_factor: f64) -> f64 {
    if mach < 0.8 {
        // Subsonic - minimal correction needed
        1.0
    } else if mach < 1.2 {
        // Transonic region - significant instability
        // Smooth transition using cosine interpolation
        let transonic_ratio = (mach - 0.8) / 0.4;
        let base_correction = 1.0 + (transonic_factor - 1.0) * (1.0 - (transonic_ratio * PI).cos()) / 2.0;
        base_correction
    } else {
        // Supersonic - stable again but reduced effect
        0.85 + 0.15 * (2.5 - mach).max(0.0) / 1.3
    }
}

/// Yaw damping factor based on stability and time
fn calculate_yaw_damping_factor(stability: f64, time_s: f64, damping_coeff: f64) -> f64 {
    // Yaw oscillations damp out over time
    // Higher stability = faster damping
    let damping_rate = damping_coeff * stability.sqrt();
    let damped = 1.0 - (-damping_rate * time_s).exp();
    
    // Ensure reasonable bounds
    damped.max(0.5).min(1.0)
}

/// Calculate equilibrium yaw angle using advanced model
pub fn calculate_advanced_yaw_of_repose(
    stability_factor: f64,
    velocity_mps: f64,
    crosswind_mps: f64,
    spin_rate_rad_s: f64,
    air_density_kg_m3: f64,
    caliber_m: f64,
) -> f64 {
    if stability_factor <= 1.0 || velocity_mps <= 0.0 {
        return 0.0;
    }
    
    // Base yaw from crosswind
    let wind_yaw = if crosswind_mps != 0.0 && velocity_mps > 0.0 {
        (crosswind_mps / velocity_mps).atan()
    } else {
        // Natural yaw from trajectory curvature (gravity-induced)
        // Empirical value based on typical trajectories
        0.001 + 0.0005 * (velocity_mps / 800.0).min(2.0)
    };
    
    // Stability-based damping (McCoy's model)
    let stability_term = ((stability_factor - 1.0) / stability_factor).sqrt();
    
    // Dynamic pressure effect
    let q = 0.5 * air_density_kg_m3 * velocity_mps.powi(2);
    let q_factor = (q / 50000.0).min(1.5).max(0.5); // Normalize around typical q
    
    // Spin effect on yaw response
    let spin_factor = if spin_rate_rad_s > 0.0 {
        let spin_param = spin_rate_rad_s * caliber_m / (2.0 * velocity_mps);
        1.0 + 0.2 * spin_param.min(0.5)
    } else {
        1.0
    };
    
    wind_yaw * stability_term * q_factor * spin_factor
}

/// Data-driven correction factor (placeholder for ML integration)
pub fn apply_ml_correction(
    base_drift: f64,
    stability: f64,
    mach: f64,
    time_s: f64,
    caliber_inches: f64,
    mass_grains: f64,
) -> f64 {
    // This function would integrate with ML models trained on real-world data
    // For now, returns the base drift unmodified
    // 
    // In production, this would:
    // 1. Extract features: [stability, mach, time_s, caliber_inches, mass_grains]
    // 2. Pass to trained neural network or gradient boosting model
    // 3. Return correction factor (typically 0.8-1.2)
    // 4. Multiply base_drift by correction factor
    
    // Placeholder implementation with simple heuristics
    let mut correction = 1.0;
    
    // Known adjustments from field data
    if stability > 2.5 && mach < 1.0 {
        correction *= 0.92; // Over-stabilized subsonic tends to drift less
    }
    
    if time_s > 2.0 && mach < 0.9 {
        correction *= 1.08; // Long flight subsonic needs more correction
    }
    
    if caliber_inches < 0.264 && mass_grains < 100.0 {
        correction *= 0.88; // Light, small caliber bullets drift less
    }
    
    base_drift * correction
}

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

    #[test]
    fn test_advanced_spin_drift() {
        // Test with typical .308 Match bullet
        let drift = calculate_advanced_spin_drift(
            1.5,        // stability
            1.2,        // time of flight
            600.0,      // current velocity m/s
            850.0,      // muzzle velocity m/s
            1500.0,     // spin rate rad/s
            0.00308,    // caliber in meters
            0.0108,     // mass in kg (168 grains)
            1.225,      // air density
            true,       // right twist
            "match",    // bullet type
        );
        
        // Should give reasonable drift (2-8 inches at 1000 yards typical)
        assert!(drift > 0.0);
        assert!(drift < 0.3); // Less than 12 inches in meters
    }

    #[test]
    fn test_transonic_correction() {
        let subsonic = calculate_transonic_correction(0.7, 0.75);
        let transonic = calculate_transonic_correction(1.0, 0.75);
        let supersonic = calculate_transonic_correction(1.5, 0.75);
        
        assert_eq!(subsonic, 1.0);
        assert!(transonic > 0.8 && transonic < 1.0);
        assert!(supersonic > 0.7 && supersonic < 1.0);
    }
}