ballistics-engine 0.14.1

use crate::transonic_drag::{get_projectile_shape, transonic_correction, ProjectileShape};
use crate::DragModel;
use ndarray::ArrayD;
use std::sync::LazyLock;
/// Drag coefficient calculations for ballistics using actual drag table data
use std::path::Path;

/// Drag table data structure
#[derive(Debug, Clone)]
pub struct DragTable {
    pub mach_values: Vec<f64>,
    pub cd_values: Vec<f64>,
}

impl DragTable {
    /// Create a new drag table from mach and cd arrays
    pub fn new(mach_values: Vec<f64>, cd_values: Vec<f64>) -> Self {
        Self {
            mach_values,
            cd_values,
        }
    }

    /// Interpolate drag coefficient for given Mach number using cubic spline-like interpolation
    pub fn interpolate(&self, mach: f64) -> f64 {
        let n = self.mach_values.len();

        if n == 0 {
            return 0.5; // Fallback
        }

        if n == 1 {
            return self.cd_values[0];
        }

        // Handle out-of-bounds cases with extrapolation
        if mach <= self.mach_values[0] {
            // Linear extrapolation below range, but prevent negative values
            if n >= 2 {
                let slope = (self.cd_values[1] - self.cd_values[0])
                    / (self.mach_values[1] - self.mach_values[0]);
                let extrapolated = self.cd_values[0] + slope * (mach - self.mach_values[0]);
                // Clamp to minimum reasonable value to prevent negative drag coefficients
                return extrapolated.max(0.01);
            }
            return self.cd_values[0];
        }

        if mach >= self.mach_values[n - 1] {
            // Linear extrapolation above range, but prevent negative values
            if n >= 2 {
                let slope = (self.cd_values[n - 1] - self.cd_values[n - 2])
                    / (self.mach_values[n - 1] - self.mach_values[n - 2]);
                let extrapolated = self.cd_values[n - 1] + slope * (mach - self.mach_values[n - 1]);
                // Clamp to minimum reasonable value to prevent negative drag coefficients
                return extrapolated.max(0.01);
            }
            return self.cd_values[n - 1];
        }

        // Find the segment containing the mach value
        let mut idx = 0;
        for i in 0..n - 1 {
            if mach >= self.mach_values[i] && mach <= self.mach_values[i + 1] {
                idx = i;
                break;
            }
        }

        // Use cubic interpolation if we have enough points, otherwise linear
        if idx > 0 && idx < n - 2 {
            // Cubic interpolation using 4 points
            self.cubic_interpolate(mach, idx)
        } else {
            // Linear interpolation for edge cases
            self.linear_interpolate(mach, idx)
        }
    }

    /// Linear interpolation between two points
    pub fn linear_interpolate(&self, mach: f64, idx: usize) -> f64 {
        // Bounds check
        if idx + 1 >= self.mach_values.len() || idx + 1 >= self.cd_values.len() {
            return self.cd_values.get(idx).copied().unwrap_or(0.5);
        }

        let x0 = self.mach_values[idx];
        let x1 = self.mach_values[idx + 1];
        let y0 = self.cd_values[idx];
        let y1 = self.cd_values[idx + 1];

        if (x1 - x0).abs() < crate::constants::MIN_DIVISION_THRESHOLD {
            return y0;
        }

        let t = (mach - x0) / (x1 - x0);
        y0 + t * (y1 - y0)
    }

    /// Cubic interpolation using 4 points (similar to cubic spline)
    pub fn cubic_interpolate(&self, mach: f64, idx: usize) -> f64 {
        // Ensure we have enough points for cubic interpolation
        if idx == 0 || idx + 1 >= self.mach_values.len() || idx + 1 >= self.cd_values.len() {
            // Fall back to linear interpolation if not enough points
            return self.linear_interpolate(mach, idx);
        }

        // Use points at idx-1, idx, idx+1, idx+2
        let x = [
            self.mach_values[idx - 1],
            self.mach_values[idx],
            self.mach_values[idx + 1],
            if idx + 2 < self.mach_values.len() {
                self.mach_values[idx + 2]
            } else {
                self.mach_values[idx + 1]
            },
        ];
        let y = [
            self.cd_values[idx - 1],
            self.cd_values[idx],
            self.cd_values[idx + 1],
            if idx + 2 < self.cd_values.len() {
                self.cd_values[idx + 2]
            } else {
                self.cd_values[idx + 1]
            },
        ];

        // Catmull-Rom spline interpolation
        // Ensure denominator is not zero
        let denominator = x[2] - x[1];
        if denominator.abs() < crate::constants::MIN_DIVISION_THRESHOLD {
            return y[1]; // Return the value at the current point if points are too close
        }
        let t = (mach - x[1]) / denominator;
        let t2 = t * t;
        let t3 = t2 * t;

        let a0 = -0.5 * y[0] + 1.5 * y[1] - 1.5 * y[2] + 0.5 * y[3];
        let a1 = y[0] - 2.5 * y[1] + 2.0 * y[2] - 0.5 * y[3];
        let a2 = -0.5 * y[0] + 0.5 * y[2];
        let a3 = y[1];

        a0 * t3 + a1 * t2 + a2 * t + a3
    }
}

/// Load drag table from NumPy binary file or CSV fallback
pub fn load_drag_table(
    drag_tables_dir: &Path,
    filename: &str,
    fallback_data: &[(f64, f64)],
) -> DragTable {
    // Try to load NumPy binary file first
    let npy_path = drag_tables_dir.join(format!("{filename}.npy"));
    if let Ok(array) = ndarray_npy::read_npy::<_, ArrayD<f64>>(&npy_path) {
        if let Ok(array_2d) = array.into_dimensionality::<ndarray::Ix2>() {
            let mach_values: Vec<f64> = array_2d.column(0).to_vec();
            let cd_values: Vec<f64> = array_2d.column(1).to_vec();
            return DragTable::new(mach_values, cd_values);
        }
    }

    // Fallback to CSV file
    let csv_path = drag_tables_dir.join(format!("{filename}.csv"));
    if let Ok(mut reader) = csv::Reader::from_path(&csv_path) {
        let mut mach_values = Vec::new();
        let mut cd_values = Vec::new();

        for record in reader.records().flatten() {
            if record.len() >= 2 {
                if let (Ok(mach), Ok(cd)) = (record[0].parse::<f64>(), record[1].parse::<f64>())
                {
                    mach_values.push(mach);
                    cd_values.push(cd);
                }
            }
        }

        if !mach_values.is_empty() {
            return DragTable::new(mach_values, cd_values);
        }
    }

    // Use fallback data if both file loading methods fail
    let mach_values: Vec<f64> = fallback_data.iter().map(|(m, _)| *m).collect();
    let cd_values: Vec<f64> = fallback_data.iter().map(|(_, cd)| *cd).collect();
    DragTable::new(mach_values, cd_values)
}

/// Find the drag tables directory relative to the current location
fn find_drag_tables_dir() -> Option<std::path::PathBuf> {
    // Try common relative paths from the Rust crate location
    let candidates = [
        "../drag_tables",
        "../../drag_tables",
        "../../../drag_tables",
        "drag_tables",
    ];

    for candidate in &candidates {
        let path = Path::new(candidate);
        if path.exists() && path.is_dir() {
            return Some(path.to_path_buf());
        }
    }

    None
}

/// G1 drag table with lazy loading
static G1_DRAG_TABLE: LazyLock<DragTable> = LazyLock::new(|| {
    let fallback_data = [
        (0.0, 0.2629),
        (0.5, 0.2695),
        (0.6, 0.2752),
        (0.7, 0.2817),
        (0.8, 0.2902),
        (0.9, 0.3012),
        (1.0, 0.4805),
        (1.1, 0.5933),
        (1.2, 0.6318),
        (1.3, 0.6440),
        (1.4, 0.6444),
        (1.5, 0.6372),
        (1.6, 0.6252),
        (1.7, 0.6105),
        (1.8, 0.5956),
        (1.9, 0.5815),
        (2.0, 0.5934),
        (2.5, 0.5598),
        (3.0, 0.5133),
        (4.0, 0.4811),
        (5.0, 0.4988),
    ];

    if let Some(drag_dir) = find_drag_tables_dir() {
        load_drag_table(&drag_dir, "g1", &fallback_data)
    } else {
        // Use fallback data if directory not found
        let mach_values: Vec<f64> = fallback_data.iter().map(|(m, _)| *m).collect();
        let cd_values: Vec<f64> = fallback_data.iter().map(|(_, cd)| *cd).collect();
        DragTable::new(mach_values, cd_values)
    }
});

/// G7 drag table with lazy loading
static G7_DRAG_TABLE: LazyLock<DragTable> = LazyLock::new(|| {
    let fallback_data = [
        (0.0, 0.1198),
        (0.5, 0.1197),
        (0.6, 0.1202),
        (0.7, 0.1213),
        (0.8, 0.1240),
        (0.9, 0.1294),
        (1.0, 0.3803),
        (1.1, 0.4015),
        (1.2, 0.4043),
        (1.3, 0.3956),
        (1.4, 0.3814),
        (1.5, 0.3663),
        (1.6, 0.3520),
        (1.7, 0.3398),
        (1.8, 0.3297),
        (1.9, 0.3221),
        (2.0, 0.2980),
        (2.5, 0.2731),
        (3.0, 0.2424),
        (4.0, 0.2196),
        (5.0, 0.1618),
    ];

    if let Some(drag_dir) = find_drag_tables_dir() {
        load_drag_table(&drag_dir, "g7", &fallback_data)
    } else {
        // Use fallback data if directory not found
        let mach_values: Vec<f64> = fallback_data.iter().map(|(m, _)| *m).collect();
        let cd_values: Vec<f64> = fallback_data.iter().map(|(_, cd)| *cd).collect();
        DragTable::new(mach_values, cd_values)
    }
});

/// G6 drag table - flat-base with 6 caliber secant ogive (military FMJ bullets)
/// MBA-156: Added for completeness with ballistics_rust
static G6_DRAG_TABLE: LazyLock<DragTable> = LazyLock::new(|| {
    let fallback_data = [
        (0.0, 0.2617),
        (0.05, 0.2553),
        (0.10, 0.2491),
        (0.15, 0.2432),
        (0.20, 0.2376),
        (0.25, 0.2324),
        (0.30, 0.2278),
        (0.35, 0.2238),
        (0.40, 0.2205),
        (0.45, 0.2177),
        (0.50, 0.2155),
        (0.55, 0.2138),
        (0.60, 0.2126),
        (0.65, 0.2121),
        (0.70, 0.2122),
        (0.75, 0.2132),
        (0.80, 0.2154),
        (0.85, 0.2194),
        (0.875, 0.2229),
        (0.90, 0.2297),
        (0.925, 0.2449),
        (0.95, 0.2732),
        (0.975, 0.3141),
        (1.0, 0.3597),
        (1.025, 0.3994),
        (1.05, 0.4261),
        (1.075, 0.4402),
        (1.10, 0.4465),
        (1.125, 0.4490),
        (1.15, 0.4497),
        (1.175, 0.4494),
        (1.20, 0.4482),
        (1.225, 0.4464),
        (1.25, 0.4441),
        (1.30, 0.4390),
        (1.35, 0.4336),
        (1.40, 0.4279),
        (1.45, 0.4221),
        (1.50, 0.4162),
        (1.55, 0.4102),
        (1.60, 0.4042),
        (1.65, 0.3981),
        (1.70, 0.3919),
        (1.75, 0.3855),
        (1.80, 0.3788),
        (1.85, 0.3721),
        (1.90, 0.3652),
        (1.95, 0.3583),
        (2.0, 0.3515),
        (2.05, 0.3447),
        (2.10, 0.3381),
        (2.15, 0.3314),
        (2.20, 0.3249),
        (2.25, 0.3185),
        (2.30, 0.3122),
        (2.35, 0.3060),
        (2.40, 0.3000),
        (2.45, 0.2941),
        (2.50, 0.2883),
        (2.60, 0.2772),
        (2.70, 0.2668),
        (2.80, 0.2574),
        (2.90, 0.2487),
        (3.0, 0.2407),
        (3.10, 0.2333),
        (3.20, 0.2265),
        (3.30, 0.2202),
        (3.40, 0.2144),
        (3.50, 0.2089),
        (3.60, 0.2039),
        (3.70, 0.1991),
        (3.80, 0.1947),
        (3.90, 0.1905),
        (4.0, 0.1866),
        (4.20, 0.1794),
        (4.40, 0.1730),
        (4.60, 0.1673),
        (4.80, 0.1621),
        (5.0, 0.1574),
    ];

    if let Some(drag_dir) = find_drag_tables_dir() {
        load_drag_table(&drag_dir, "g6", &fallback_data)
    } else {
        // Use fallback data if directory not found
        let mach_values: Vec<f64> = fallback_data.iter().map(|(m, _)| *m).collect();
        let cd_values: Vec<f64> = fallback_data.iter().map(|(_, cd)| *cd).collect();
        DragTable::new(mach_values, cd_values)
    }
});

/// G8 drag table - flat-base with 10 caliber secant ogive
/// MBA-156: Added for completeness with ballistics_rust
static G8_DRAG_TABLE: LazyLock<DragTable> = LazyLock::new(|| {
    let fallback_data = [
        (0.0, 0.2105),
        (0.05, 0.2105),
        (0.10, 0.2104),
        (0.15, 0.2104),
        (0.20, 0.2103),
        (0.25, 0.2103),
        (0.30, 0.2103),
        (0.35, 0.2103),
        (0.40, 0.2103),
        (0.45, 0.2102),
        (0.50, 0.2102),
        (0.55, 0.2102),
        (0.60, 0.2102),
        (0.65, 0.2102),
        (0.70, 0.2103),
        (0.75, 0.2103),
        (0.80, 0.2104),
        (0.825, 0.2104),
        (0.85, 0.2105),
        (0.875, 0.2106),
        (0.90, 0.2109),
        (0.925, 0.2183),
        (0.95, 0.2571),
        (0.975, 0.3358),
        (1.0, 0.4068),
        (1.025, 0.4378),
        (1.05, 0.4476),
        (1.075, 0.4493),
        (1.10, 0.4477),
        (1.125, 0.4450),
        (1.15, 0.4419),
        (1.20, 0.4353),
        (1.25, 0.4283),
        (1.30, 0.4208),
        (1.35, 0.4133),
        (1.40, 0.4059),
        (1.45, 0.3986),
        (1.50, 0.3915),
        (1.55, 0.3845),
        (1.60, 0.3777),
        (1.65, 0.3710),
        (1.70, 0.3645),
        (1.75, 0.3581),
        (1.80, 0.3519),
        (1.85, 0.3458),
        (1.90, 0.3400),
        (1.95, 0.3343),
        (2.0, 0.3288),
        (2.05, 0.3234),
        (2.10, 0.3182),
        (2.15, 0.3131),
        (2.20, 0.3081),
        (2.25, 0.3032),
        (2.30, 0.2983),
        (2.35, 0.2937),
        (2.40, 0.2891),
        (2.45, 0.2845),
        (2.50, 0.2802),
        (2.60, 0.2720),
        (2.70, 0.2642),
        (2.80, 0.2569),
        (2.90, 0.2499),
        (3.0, 0.2432),
        (3.10, 0.2368),
        (3.20, 0.2308),
        (3.30, 0.2251),
        (3.40, 0.2197),
        (3.50, 0.2147),
        (3.60, 0.2101),
        (3.70, 0.2058),
        (3.80, 0.2019),
        (3.90, 0.1983),
        (4.0, 0.1950),
        (4.20, 0.1890),
        (4.40, 0.1837),
        (4.60, 0.1791),
        (4.80, 0.1750),
        (5.0, 0.1713),
    ];

    if let Some(drag_dir) = find_drag_tables_dir() {
        load_drag_table(&drag_dir, "g8", &fallback_data)
    } else {
        // Use fallback data if directory not found
        let mach_values: Vec<f64> = fallback_data.iter().map(|(m, _)| *m).collect();
        let cd_values: Vec<f64> = fallback_data.iter().map(|(_, cd)| *cd).collect();
        DragTable::new(mach_values, cd_values)
    }
});

/// Get drag coefficient for given Mach number and drag model
pub fn get_drag_coefficient(mach: f64, drag_model: &DragModel) -> f64 {
    match drag_model {
        DragModel::G1 => G1_DRAG_TABLE.interpolate(mach),
        DragModel::G6 => G6_DRAG_TABLE.interpolate(mach),
        DragModel::G7 => G7_DRAG_TABLE.interpolate(mach),
        DragModel::G8 => G8_DRAG_TABLE.interpolate(mach),
        _ => G1_DRAG_TABLE.interpolate(mach), // Default to G1 for other models
    }
}

/// Get drag coefficient with optional transonic corrections
pub fn get_drag_coefficient_with_transonic(
    mach: f64,
    drag_model: &DragModel,
    apply_transonic_correction: bool,
    projectile_shape: Option<ProjectileShape>,
    caliber: Option<f64>,
    weight_grains: Option<f64>,
) -> f64 {
    // Get base drag coefficient
    let base_cd = get_drag_coefficient(mach, drag_model);

    // Apply transonic corrections if requested and in transonic regime
    if apply_transonic_correction && (0.8..=1.3).contains(&mach) {
        // Determine projectile shape if not provided
        let shape = match projectile_shape {
            Some(s) => s,
            None => {
                if let (Some(cal), Some(weight)) = (caliber, weight_grains) {
                    get_projectile_shape(
                        cal,
                        weight,
                        match drag_model {
                            DragModel::G1 => "G1",
                            DragModel::G6 => "G6",
                            DragModel::G7 => "G7",
                            DragModel::G8 => "G8",
                            _ => "G1", // Default to G1
                        },
                    )
                } else {
                    ProjectileShape::Spitzer // Default
                }
            }
        };

        // Apply correction
        transonic_correction(mach, base_cd, shape, true)
    } else {
        base_cd
    }
}

/// Get drag coefficient with optional transonic and Reynolds corrections
pub fn get_drag_coefficient_full(
    mach: f64,
    drag_model: &DragModel,
    apply_transonic_correction: bool,
    apply_reynolds_correction: bool,
    projectile_shape: Option<ProjectileShape>,
    caliber: Option<f64>,
    weight_grains: Option<f64>,
    velocity_mps: Option<f64>,
    air_density_kg_m3: Option<f64>,
    temperature_c: Option<f64>,
) -> f64 {
    // Get base drag coefficient with transonic corrections if applicable
    let mut cd = get_drag_coefficient_with_transonic(
        mach,
        drag_model,
        apply_transonic_correction,
        projectile_shape,
        caliber,
        weight_grains,
    );

    // Apply Reynolds corrections for low velocities (subsonic only)
    if apply_reynolds_correction && mach < 1.0 {
        if let (Some(v), Some(cal), Some(rho), Some(temp)) =
            (velocity_mps, caliber, air_density_kg_m3, temperature_c)
        {
            use crate::reynolds::apply_reynolds_correction;
            cd = apply_reynolds_correction(cd, v, cal, rho, temp, mach);
        }
    }

    cd
}

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

    #[test]
    fn test_g1_drag_coefficient_interpolation() {
        let cd = get_drag_coefficient(1.0, &DragModel::G1);
        // Should be close to the G1 standard value at Mach 1.0
        assert!(cd > 0.4 && cd < 0.6, "G1 CD at Mach 1.0: {cd}");
    }

    #[test]
    fn test_g7_drag_coefficient_interpolation() {
        let cd = get_drag_coefficient(1.0, &DragModel::G7);
        // Should be close to the G7 standard value at Mach 1.0
        assert!(cd > 0.3 && cd < 0.5, "G7 CD at Mach 1.0: {cd}");
    }

    #[test]
    fn test_drag_coefficient_continuity() {
        // Test that drag coefficient function is smooth
        for mach in [0.5, 0.8, 1.0, 1.2, 1.5, 2.0, 3.0] {
            let cd_before = get_drag_coefficient(mach - 0.01, &DragModel::G1);
            let cd_after = get_drag_coefficient(mach + 0.01, &DragModel::G1);
            let difference = (cd_after - cd_before).abs();
            assert!(
                difference < 0.05,
                "Large discontinuity at Mach {mach}: {cd_before} vs {cd_after}"
            );
        }
    }

    #[test]
    fn test_extrapolation_bounds() {
        // Test extrapolation below range
        let cd_low = get_drag_coefficient(0.0, &DragModel::G1);
        assert!(cd_low > 0.01 && cd_low < 0.5, "Low Mach G1: {cd_low}");

        // Test extrapolation above range - should be clamped to positive values
        let cd_high = get_drag_coefficient(10.0, &DragModel::G1);
        assert!(cd_high > 0.01, "High Mach G1 should be positive: {cd_high}");

        // Same for G7
        let cd_low_g7 = get_drag_coefficient(0.0, &DragModel::G7);
        assert!(
            cd_low_g7 > 0.01,
            "Low Mach G7 should be positive: {cd_low_g7}"
        );

        let cd_high_g7 = get_drag_coefficient(20.0, &DragModel::G7);
        assert!(
            cd_high_g7 >= 0.01,
            "High Mach G7 should be positive: {cd_high_g7}"
        );
    }

    #[test]
    fn test_drag_table_creation() {
        let mach_vals = vec![0.5, 1.0, 1.5, 2.0];
        let cd_vals = vec![0.2, 0.5, 0.4, 0.3];
        let table = DragTable::new(mach_vals, cd_vals);

        // Test exact interpolation
        assert!((table.interpolate(1.0) - 0.5).abs() < 1e-10);

        // Test interpolation between points
        let cd_interp = table.interpolate(1.25);
        assert!(cd_interp > 0.4 && cd_interp < 0.5);
    }

    #[test]
    fn test_drag_table_empty() {
        let table = DragTable::new(vec![], vec![]);
        let result = table.interpolate(1.0);
        assert_eq!(result, 0.5); // Should return fallback value
    }

    #[test]
    fn test_drag_table_single_point() {
        let table = DragTable::new(vec![1.0], vec![0.4]);

        // Should return the single value for any Mach
        assert_eq!(table.interpolate(0.5), 0.4);
        assert_eq!(table.interpolate(1.0), 0.4);
        assert_eq!(table.interpolate(2.0), 0.4);
    }

    #[test]
    fn test_drag_table_two_points() {
        let table = DragTable::new(vec![1.0, 2.0], vec![0.4, 0.6]);

        // Exact matches
        assert!((table.interpolate(1.0) - 0.4).abs() < 1e-10);
        assert!((table.interpolate(2.0) - 0.6).abs() < 1e-10);

        // Linear interpolation
        let mid = table.interpolate(1.5);
        assert!((mid - 0.5).abs() < 1e-10);

        // Extrapolation
        let below = table.interpolate(0.5);
        assert!(below.abs() > 1e-10); // Should have some value

        let above = table.interpolate(3.0);
        assert!(above.abs() > 1e-10); // Should have some value
    }

    #[test]
    fn test_linear_interpolation() {
        let table = DragTable::new(vec![0.0, 1.0, 2.0], vec![0.2, 0.5, 0.3]);

        // Test linear interpolation between first two points
        let result = table.linear_interpolate(0.5, 0);
        assert!((result - 0.35).abs() < 1e-10);

        // Test edge case with zero denominator
        let table_same = DragTable::new(vec![1.0, 1.0], vec![0.4, 0.6]);
        let result_same = table_same.linear_interpolate(1.0, 0);
        assert_eq!(result_same, 0.4); // Should return first value
    }

    #[test]
    fn test_cubic_interpolation() {
        // Create a table with enough points for cubic interpolation
        let table = DragTable::new(vec![0.5, 1.0, 1.5, 2.0, 2.5], vec![0.2, 0.4, 0.6, 0.5, 0.3]);

        // Test cubic interpolation in the middle
        let result = table.cubic_interpolate(1.25, 1);

        // Should be between the neighboring values
        assert!(result > 0.3 && result < 0.7);

        // Should be smooth (not exactly linear)
        let linear_result = table.linear_interpolate(1.25, 1);
        // Cubic and linear should be close but not identical for smooth curves
        assert!((result - linear_result).abs() < 0.2);
    }

    #[test]
    fn test_find_drag_tables_dir() {
        // This test may pass or fail depending on the environment
        // but should not panic
        let _dir = find_drag_tables_dir();
        // Just ensure the function doesn't panic
    }

    #[test]
    fn test_load_drag_table_fallback() {
        use std::path::Path;

        // Test with non-existent directory - should use fallback data
        let fake_dir = Path::new("/non/existent/directory");
        let fallback_data = [(0.5, 0.2), (1.0, 0.4), (1.5, 0.3)];

        let table = load_drag_table(fake_dir, "test", &fallback_data);

        // Should have fallback data
        assert_eq!(table.mach_values.len(), 3);
        assert_eq!(table.cd_values.len(), 3);
        assert_eq!(table.mach_values[0], 0.5);
        assert_eq!(table.cd_values[0], 0.2);
    }

    #[test]
    fn test_known_drag_values() {
        // Test against known ballistic standard values

        // G1 at Mach 1.0 should be around 0.4805
        let g1_mach1 = get_drag_coefficient(1.0, &DragModel::G1);
        assert!(
            (g1_mach1 - 0.4805).abs() < 0.01,
            "G1 at Mach 1.0: {g1_mach1}"
        );

        // G7 at Mach 1.0 should be around 0.3803
        let g7_mach1 = get_drag_coefficient(1.0, &DragModel::G7);
        assert!(
            (g7_mach1 - 0.3803).abs() < 0.01,
            "G7 at Mach 1.0: {g7_mach1}"
        );

        // G1 should generally be higher than G7 in transonic region
        assert!(g1_mach1 > g7_mach1, "G1 should be > G7 at Mach 1.0");
    }

    #[test]
    fn test_monotonicity_properties() {
        // Test general drag curve properties

        // G1 should peak somewhere in transonic region
        let mach_values: Vec<f64> = (8..20).map(|i| i as f64 * 0.1).collect(); // 0.8 to 1.9
        let g1_values: Vec<f64> = mach_values
            .iter()
            .map(|&m| get_drag_coefficient(m, &DragModel::G1))
            .collect();

        // Find maximum
        let max_value = g1_values.iter().copied().fold(0.0_f64, f64::max);
        let max_index = g1_values
            .iter()
            .position(|&x| x == max_value)
            .expect("Should find maximum in non-empty vector");
        let peak_mach = mach_values
            .get(max_index)
            .copied()
            .expect("Index should be valid");

        // Peak should be in reasonable range
        assert!(
            peak_mach > 1.0 && peak_mach < 1.6,
            "G1 peak at Mach {peak_mach}"
        );
        assert!(
            max_value > 0.5 && max_value < 1.0,
            "G1 peak value: {max_value}"
        );
    }

    #[test]
    fn test_physical_constraints() {
        let test_machs = [0.1, 0.5, 0.8, 1.0, 1.2, 1.5, 2.0, 3.0, 5.0];

        for &mach in &test_machs {
            let g1_cd = get_drag_coefficient(mach, &DragModel::G1);
            let g7_cd = get_drag_coefficient(mach, &DragModel::G7);

            // All drag coefficients should be positive
            assert!(g1_cd > 0.0, "G1 CD negative at Mach {mach}: {g1_cd}");
            assert!(g7_cd > 0.0, "G7 CD negative at Mach {mach}: {g7_cd}");

            // Should be in reasonable physical ranges
            assert!(g1_cd < 2.0, "G1 CD too high at Mach {mach}: {g1_cd}");
            assert!(g7_cd < 1.5, "G7 CD too high at Mach {mach}: {g7_cd}");
        }
    }

    #[test]
    fn test_performance_characteristics() {
        // This test ensures the implementation is efficient
        use std::time::Instant;

        let start = Instant::now();

        // Perform many calculations
        for i in 0..1000 {
            let mach = 0.5 + (i as f64) * 0.004; // 0.5 to 4.5
            let _g1 = get_drag_coefficient(mach, &DragModel::G1);
            let _g7 = get_drag_coefficient(mach, &DragModel::G7);
        }

        let elapsed = start.elapsed();

        // Should complete 2000 calculations quickly (within 100ms)
        assert!(
            elapsed.as_millis() < 100,
            "Performance test took {}ms",
            elapsed.as_millis()
        );
    }
}

/// Interpolate BC value for given Mach number from segments
pub fn interpolated_bc(mach: f64, segments: &[(f64, f64)]) -> f64 {
    if segments.is_empty() {
        return crate::constants::BC_FALLBACK_CONSERVATIVE; // Conservative fallback based on database analysis
    }

    // Get just the mach values
    let mach_values: Vec<f64> = segments.iter().map(|(m, _)| *m).collect();

    // Double-check we have values after collection
    if mach_values.is_empty() || segments.is_empty() {
        return crate::constants::BC_FALLBACK_CONSERVATIVE; // Conservative fallback based on database analysis
    }

    // Handle edge cases with safe indexing
    if let Some(first_mach) = mach_values.first() {
        if mach <= *first_mach {
            return segments.first().map(|(_, bc)| *bc).unwrap_or(0.5);
        }
    }

    if let Some(last_mach) = mach_values.last() {
        if mach >= *last_mach {
            return segments.last().map(|(_, bc)| *bc).unwrap_or(0.5);
        }
    }

    // Binary search to find the right segment with safe comparison
    let idx = match mach_values
        .binary_search_by(|&m| m.partial_cmp(&mach).unwrap_or(std::cmp::Ordering::Equal))
    {
        Ok(idx) => {
            // Exact match - safely get the BC value
            return segments.get(idx).map(|(_, bc)| *bc).unwrap_or(0.5);
        }
        Err(idx) => idx, // Insert position
    };

    // Ensure idx is valid for interpolation
    if idx == 0 || idx >= segments.len() {
        // Shouldn't happen given the edge case checks above, but be defensive
        // Use safe indexing
        let safe_idx = idx.saturating_sub(1).min(segments.len().saturating_sub(1));
        return segments.get(safe_idx).map(|(_, bc)| *bc).unwrap_or(0.5);
    }

    // Linear interpolation between the two closest points with safe indexing
    match (segments.get(idx - 1), segments.get(idx)) {
        (Some((lo_mach, lo_bc)), Some((hi_mach, hi_bc))) => {
            // Ensure denominator is not zero for safe interpolation
            let denominator = hi_mach - lo_mach;
            if denominator.abs() < crate::constants::MIN_DIVISION_THRESHOLD {
                return *lo_bc; // Return lower BC if Mach values are too close
            }
            let frac = (mach - lo_mach) / denominator;
            lo_bc + frac * (hi_bc - lo_bc)
        }
        _ => 0.5, // Fallback if indices are somehow invalid
    }
}

// Removed Python-specific function