ballistics_engine/

bc_table_5d.rs

// BC5D - 5-Dimensional BC Correction Table with Caliber-Specific Files
//
// This module provides offline BC corrections by loading precomputed tables
// of correction factors derived from ML model predictions. The tables are
// caliber-specific and indexed by:
//   - Weight (grains) - caliber-specific ranges
//   - Base BC (0.05-1.2)
//   - Muzzle Velocity (2000-4000 fps)
//   - Current Velocity (500-4000 fps, dense in transonic)
//   - Drag Model (G1, G7)
//
// Binary file format (BC5D v2):
//   Header (80 bytes):
//     - Magic: 4 bytes ('BC5D')
//     - Version: 4 bytes (uint32)
//     - Caliber: 4 bytes (float32)
//     - Flags: 4 bytes (uint32)
//     - Padding: 4 bytes
//     - dim_weight: 4 bytes (uint32)
//     - dim_bc: 4 bytes (uint32)
//     - dim_muzzle_vel: 4 bytes (uint32)
//     - dim_current_vel: 4 bytes (uint32)
//     - dim_drag_types: 4 bytes (uint32)
//     - timestamp: 8 bytes (uint64)
//     - checksum: 4 bytes (uint32, CRC32 of data section)
//     - api_version: 16 bytes (null-padded string)
//     - reserved: 12 bytes
//   Bin definitions:
//     - Weight bins: dim_weight * 4 bytes (float32)
//     - BC bins: dim_bc * 4 bytes (float32)
//     - Muzzle velocity bins: dim_muzzle_vel * 4 bytes (float32)
//     - Current velocity bins: dim_current_vel * 4 bytes (float32)
//   Data section:
//     - Correction factors: total_cells * 4 bytes (float32)
//     - Layout: [drag_type][weight][bc][muzzle_vel][current_vel]
//
// Correction factors are ratios: predicted_bc / base_bc
// Range: 0.5 to 1.5 (clipped during generation)

use std::collections::HashMap;
use std::fs::File;
use std::io::{BufReader, Read};
use std::path::{Path, PathBuf};

/// Magic bytes for BC5D format
const MAGIC: &[u8; 4] = b"BC5D";

/// Supported format version
const SUPPORTED_VERSION: u32 = 2;

/// BC5D table with 4D interpolation (drag type is discrete)
#[derive(Debug)]
pub struct Bc5dTable {
    /// Caliber this table is for
    caliber: f32,
    /// Correction data: [drag_type][weight][bc][muzzle_vel][current_vel]
    data: Vec<f32>,
    /// Weight bin values (grains)
    weight_bins: Vec<f32>,
    /// BC bin values
    bc_bins: Vec<f32>,
    /// Muzzle velocity bin values (fps)
    muzzle_vel_bins: Vec<f32>,
    /// Current velocity bin values (fps)
    current_vel_bins: Vec<f32>,
    /// Number of drag types (typically 2: G1=0, G7=1)
    num_drag_types: usize,
    /// Table version
    version: u32,
    /// API version used to generate the table
    api_version: String,
    /// Generation timestamp
    timestamp: u64,
}

/// A velocity-keyed BC schedule and the scalar BC used for any interior coverage gap.
#[cfg(any(test, target_arch = "wasm32"))]
pub(crate) struct Bc5dSegmentSchedule {
    pub(crate) segments: Vec<crate::BCSegmentData>,
    pub(crate) fallback_bc: f64,
}

/// Manager for loading caliber-specific BC5D tables
#[derive(Debug, Default)]
pub struct Bc5dTableManager {
    /// Directory containing BC5D table files
    table_dir: Option<PathBuf>,
    /// Loaded tables by caliber (rounded to 3 decimal places)
    tables: HashMap<i32, Bc5dTable>,
}

/// Error type for BC5D table operations
#[derive(Debug)]
pub enum Bc5dError {
    IoError(std::io::Error),
    InvalidMagic,
    UnsupportedVersion(u32),
    ChecksumMismatch { expected: u32, actual: u32 },
    InvalidDimensions,
    TableNotFound(f64),
    NoTableDirectory,
}

impl std::fmt::Display for Bc5dError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            Bc5dError::IoError(e) => write!(f, "IO error: {}", e),
            Bc5dError::InvalidMagic => write!(f, "Invalid file magic (expected 'BC5D')"),
            Bc5dError::UnsupportedVersion(v) => write!(f, "Unsupported table version: {}", v),
            Bc5dError::ChecksumMismatch { expected, actual } => {
                write!(f, "Checksum mismatch: expected {:08x}, got {:08x}", expected, actual)
            }
            Bc5dError::InvalidDimensions => write!(f, "Invalid table dimensions"),
            Bc5dError::TableNotFound(cal) => write!(f, "No BC5D table found for caliber {:.3}", cal),
            Bc5dError::NoTableDirectory => write!(f, "No BC table directory configured"),
        }
    }
}

impl std::error::Error for Bc5dError {}

impl From<std::io::Error> for Bc5dError {
    fn from(e: std::io::Error) -> Self {
        Bc5dError::IoError(e)
    }
}

impl Bc5dTable {
    /// Load a BC5D table from a binary file
    pub fn load<P: AsRef<Path>>(path: P) -> Result<Self, Bc5dError> {
        let file = File::open(&path)?;
        Self::from_reader(BufReader::new(file))
    }

    /// Parse a BC5D table directly from an in-memory byte slice.
    ///
    /// Behaves identically to [`Bc5dTable::load`] but performs no filesystem
    /// access, which makes it usable from WASM (`wasm32-unknown-unknown`, where
    /// there is no `std::fs`). The host JS/Node layer fetches the `.bin` and
    /// hands the raw bytes across the boundary.
    pub fn from_bytes(bytes: &[u8]) -> Result<Self, Bc5dError> {
        Self::from_reader(std::io::Cursor::new(bytes))
    }

    fn from_reader<R: Read>(mut reader: R) -> Result<Self, Bc5dError> {
        // Read and validate magic
        let mut magic = [0u8; 4];
        reader.read_exact(&mut magic)?;
        if &magic != MAGIC {
            return Err(Bc5dError::InvalidMagic);
        }

        // Read header fields
        let version = read_u32(&mut reader)?;
        if version != SUPPORTED_VERSION {
            return Err(Bc5dError::UnsupportedVersion(version));
        }

        let caliber = read_f32(&mut reader)?;
        let _flags = read_u32(&mut reader)?;
        let _padding = read_u32(&mut reader)?;

        let dim_weight = read_u32(&mut reader)? as usize;
        let dim_bc = read_u32(&mut reader)? as usize;
        let dim_muzzle_vel = read_u32(&mut reader)? as usize;
        let dim_current_vel = read_u32(&mut reader)? as usize;
        let dim_drag_types = read_u32(&mut reader)? as usize;

        let timestamp = read_u64(&mut reader)?;
        let stored_checksum = read_u32(&mut reader)?;

        // Read API version (16 bytes, null-terminated)
        let mut api_version_bytes = [0u8; 16];
        reader.read_exact(&mut api_version_bytes)?;
        let api_version = String::from_utf8_lossy(&api_version_bytes)
            .trim_end_matches('\0')
            .to_string();

        // Skip reserved bytes
        let mut reserved = [0u8; 12];
        reader.read_exact(&mut reserved)?;

        // Validate dimensions
        if dim_weight == 0 || dim_bc == 0 || dim_muzzle_vel == 0 || dim_current_vel == 0 || dim_drag_types == 0 {
            return Err(Bc5dError::InvalidDimensions);
        }

        // Read bin definitions
        let weight_bins = read_f32_array(&mut reader, dim_weight)?;
        let bc_bins = read_f32_array(&mut reader, dim_bc)?;
        let muzzle_vel_bins = read_f32_array(&mut reader, dim_muzzle_vel)?;
        let current_vel_bins = read_f32_array(&mut reader, dim_current_vel)?;

        // Read data section. Bound the product with checked arithmetic so a corrupt or
        // hostile file cannot overflow (debug panic / release wrap) or trigger a huge OOM
        // allocation before the trailing CRC check can reject it.
        const MAX_TOTAL_CELLS: usize = 64_000_000; // ~256 MB of f32; far above any real table
        let total_cells = dim_drag_types
            .checked_mul(dim_weight)
            .and_then(|x| x.checked_mul(dim_bc))
            .and_then(|x| x.checked_mul(dim_muzzle_vel))
            .and_then(|x| x.checked_mul(dim_current_vel))
            .filter(|&n| n <= MAX_TOTAL_CELLS)
            .ok_or(Bc5dError::InvalidDimensions)?;
        let data = read_f32_array(&mut reader, total_cells)?;

        // Verify checksum (CRC32 of bins + data)
        let mut checksum_data = Vec::new();
        for &v in &weight_bins {
            checksum_data.extend_from_slice(&v.to_le_bytes());
        }
        for &v in &bc_bins {
            checksum_data.extend_from_slice(&v.to_le_bytes());
        }
        for &v in &muzzle_vel_bins {
            checksum_data.extend_from_slice(&v.to_le_bytes());
        }
        for &v in &current_vel_bins {
            checksum_data.extend_from_slice(&v.to_le_bytes());
        }
        for &v in &data {
            checksum_data.extend_from_slice(&v.to_le_bytes());
        }

        let calculated_checksum = crc32_ieee(&checksum_data);
        if calculated_checksum != stored_checksum {
            return Err(Bc5dError::ChecksumMismatch {
                expected: stored_checksum,
                actual: calculated_checksum,
            });
        }

        Ok(Bc5dTable {
            caliber,
            data,
            weight_bins,
            bc_bins,
            muzzle_vel_bins,
            current_vel_bins,
            num_drag_types: dim_drag_types,
            version,
            api_version,
            timestamp,
        })
    }

    /// Look up a BC correction factor with 4D linear interpolation
    /// (drag type is discrete, not interpolated)
    ///
    /// # Arguments
    /// * `weight_grains` - Bullet weight in grains
    /// * `base_bc` - Published BC value
    /// * `muzzle_velocity` - Initial muzzle velocity in fps
    /// * `current_velocity` - Current bullet velocity in fps
    /// * `drag_type` - "G1" or "G7"
    ///
    /// # Returns
    /// Correction factor (multiply published BC by this value)
    pub fn lookup(
        &self,
        weight_grains: f64,
        base_bc: f64,
        muzzle_velocity: f64,
        current_velocity: f64,
        drag_type: &str,
    ) -> f64 {
        // Get drag type index (0 = G1, 1 = G7)
        let drag_idx = if drag_type.eq_ignore_ascii_case("G7") { 1 } else { 0 };

        // Clamp drag_idx to valid range
        let drag_idx = drag_idx.min(self.num_drag_types - 1);

        // Find interpolation indices and weights for each continuous dimension
        let (weight_idx, weight_w) = self.interp_idx(weight_grains as f32, &self.weight_bins);
        let (bc_idx, bc_w) = self.interp_idx(base_bc as f32, &self.bc_bins);
        let (muzzle_idx, muzzle_w) = self.interp_idx(muzzle_velocity as f32, &self.muzzle_vel_bins);
        let (current_idx, current_w) = self.interp_idx(current_velocity as f32, &self.current_vel_bins);

        // 4D linear interpolation (16 corners of a hypercube)
        let mut result = 0.0f64;

        for dw in 0..2 {
            for db in 0..2 {
                for dm in 0..2 {
                    for dc in 0..2 {
                        // Calculate weight for this corner
                        let weight = (if dw == 0 { 1.0 - weight_w } else { weight_w })
                            * (if db == 0 { 1.0 - bc_w } else { bc_w })
                            * (if dm == 0 { 1.0 - muzzle_w } else { muzzle_w })
                            * (if dc == 0 { 1.0 - current_w } else { current_w });

                        // Get clamped indices
                        let wi = (weight_idx + dw).min(self.weight_bins.len() - 1);
                        let bi = (bc_idx + db).min(self.bc_bins.len() - 1);
                        let mi = (muzzle_idx + dm).min(self.muzzle_vel_bins.len() - 1);
                        let ci = (current_idx + dc).min(self.current_vel_bins.len() - 1);

                        // Calculate flat index
                        let idx = self.flat_index(drag_idx, wi, bi, mi, ci);
                        result += weight * self.data[idx] as f64;
                    }
                }
            }
        }

        // A correction is multiplicative, so an undefined table result must be neutral rather
        // than silently becoming the most aggressive allowed degradation.
        if !result.is_finite() {
            return 1.0;
        }

        result.clamp(0.5, 1.5)
    }

    /// Get the effective BC at a given velocity
    ///
    /// This multiplies the base BC by the correction factor from the table.
    pub fn get_effective_bc(
        &self,
        weight_grains: f64,
        base_bc: f64,
        muzzle_velocity: f64,
        current_velocity: f64,
        drag_type: &str,
    ) -> f64 {
        let correction = self.lookup(weight_grains, base_bc, muzzle_velocity, current_velocity, drag_type);
        base_bc * correction
    }

    /// Generate velocity-dependent BC segments for a bullet from this table.
    ///
    /// This mirrors the CLI's `--bc-table-dir` segment synthesis: the 4D
    /// correction surface is sampled at a fixed ladder of velocity breakpoints
    /// (from 500 fps up through the muzzle velocity) and each adjacent pair
    /// becomes a [`crate::BCSegmentData`] carrying the corrected BC over that
    /// band. The solver consumes these via `inputs.bc_segments_data` +
    /// `use_bc_segments`, giving the same velocity-dependent BC degradation
    /// offline that the online solver produces.
    ///
    /// All velocities are in fps and `weight_grains` in grains, matching the
    /// table's native units. Returns `None` when the table carries no
    /// meaningful correction for this bullet (every sampled cell ~= 1.0), so
    /// callers can leave the constant published BC in place.
    pub fn generate_segments(
        &self,
        base_bc: f64,
        drag_type: &str,
        weight_grains: f64,
        muzzle_velocity_fps: Option<f64>,
    ) -> Option<Vec<crate::BCSegmentData>> {
        let mut breakpoints: Vec<f64> = vec![
            4000.0, 3500.0, 3000.0, 2700.0, 2500.0, 2300.0, 2100.0, 2000.0, 1900.0, 1800.0, 1700.0,
            1600.0, 1500.0, 1400.0, 1350.0, 1300.0, 1250.0, 1200.0, 1150.0, 1100.0, 1050.0, 1000.0,
            950.0, 900.0, 850.0, 800.0, 700.0, 600.0, 500.0,
        ];
        if let Some(mv) = muzzle_velocity_fps {
            breakpoints.push(mv);
        }

        let mut velocities: Vec<f64> = breakpoints
            .into_iter()
            .filter(|&v| v >= 500.0 && muzzle_velocity_fps.map_or(true, |mv| v <= mv))
            .collect();
        velocities.sort_by(|a, b| b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal));
        velocities.dedup();

        // Correction factors were generated relative to the highest (muzzle)
        // velocity, so anchor every lookup to that reference.
        let reference_mv = velocities.first().copied().unwrap_or(3000.0);

        let mut segments: Vec<crate::BCSegmentData> = Vec::new();
        let mut any_correction = false;
        for i in 0..velocities.len().saturating_sub(1) {
            let vel_max = velocities[i];
            let vel_min = velocities[i + 1];
            let vel_mid = (vel_max + vel_min) / 2.0;

            let correction =
                self.lookup(weight_grains, base_bc, reference_mv, vel_mid, drag_type);
            if (correction - 1.0).abs() > 1e-6 {
                any_correction = true;
            }
            segments.push(crate::BCSegmentData {
                velocity_min: vel_min,
                velocity_max: vel_max,
                bc_value: base_bc * correction,
            });
        }

        if any_correction && !segments.is_empty() {
            Some(segments)
        } else {
            None
        }
    }

    /// Generate the BC5D schedule consumed by the WASM frontend.
    #[cfg(any(test, target_arch = "wasm32"))]
    pub(crate) fn generate_segment_schedule(
        &self,
        base_bc: f64,
        drag_type: &str,
        weight_grains: f64,
        muzzle_velocity_fps: f64,
    ) -> Option<Bc5dSegmentSchedule> {
        let segments =
            self.generate_segments(base_bc, drag_type, weight_grains, Some(muzzle_velocity_fps))?;
        let fallback_bc = self.get_effective_bc(
            weight_grains,
            base_bc,
            muzzle_velocity_fps,
            muzzle_velocity_fps,
            drag_type,
        );

        Some(Bc5dSegmentSchedule {
            segments,
            fallback_bc,
        })
    }

    /// Find interpolation index and weight for a value in bins
    fn interp_idx(&self, value: f32, bins: &[f32]) -> (usize, f64) {
        if bins.len() < 2 || value.is_nan() {
            return (0, 0.0);
        }

        // Handle out of range (clamp to edges)
        if value <= bins[0] {
            return (0, 0.0);
        }
        if value >= bins[bins.len() - 1] {
            return (bins.len().saturating_sub(2), 1.0);
        }

        // Binary search for interval containing value
        let last_interval = bins.len().saturating_sub(2);
        let idx = match bins.binary_search_by(|probe| {
            probe
                .partial_cmp(&value)
                .unwrap_or(std::cmp::Ordering::Equal)
        }) {
            Ok(i) => i.saturating_sub(1).min(last_interval),
            Err(i) => i.saturating_sub(1).min(last_interval),
        };

        // Calculate interpolation weight
        let low = bins[idx];
        let high = bins[idx + 1];
        let weight = if high > low {
            ((value - low) / (high - low)) as f64
        } else {
            0.0
        };

        (idx, weight)
    }

    /// Calculate flat array index from 5D indices
    fn flat_index(&self, drag_idx: usize, weight_idx: usize, bc_idx: usize, muzzle_idx: usize, current_idx: usize) -> usize {
        let n_weight = self.weight_bins.len();
        let n_bc = self.bc_bins.len();
        let n_muzzle = self.muzzle_vel_bins.len();
        let n_current = self.current_vel_bins.len();

        drag_idx * (n_weight * n_bc * n_muzzle * n_current)
            + weight_idx * (n_bc * n_muzzle * n_current)
            + bc_idx * (n_muzzle * n_current)
            + muzzle_idx * n_current
            + current_idx
    }

    /// Get caliber this table is for
    pub fn caliber(&self) -> f32 {
        self.caliber
    }

    /// Get table version
    pub fn version(&self) -> u32 {
        self.version
    }

    /// Get API version used to generate the table
    pub fn api_version(&self) -> &str {
        &self.api_version
    }

    /// Get generation timestamp
    pub fn timestamp(&self) -> u64 {
        self.timestamp
    }

    /// Get total number of cells in the table
    pub fn total_cells(&self) -> usize {
        self.data.len()
    }

    /// Get table dimensions as a string
    pub fn dimensions_str(&self) -> String {
        format!(
            "{}x{}x{}x{}x{} (weight x bc x muzzle_vel x current_vel x drag_types)",
            self.weight_bins.len(),
            self.bc_bins.len(),
            self.muzzle_vel_bins.len(),
            self.current_vel_bins.len(),
            self.num_drag_types
        )
    }

    /// Get weight range
    pub fn weight_range(&self) -> (f32, f32) {
        (*self.weight_bins.first().unwrap_or(&0.0), *self.weight_bins.last().unwrap_or(&0.0))
    }

    /// Get velocity range
    pub fn velocity_range(&self) -> (f32, f32) {
        (*self.current_vel_bins.first().unwrap_or(&0.0), *self.current_vel_bins.last().unwrap_or(&0.0))
    }
}

impl Bc5dTableManager {
    /// Create a new table manager with a directory path
    pub fn new<P: AsRef<Path>>(table_dir: P) -> Self {
        Bc5dTableManager {
            table_dir: Some(table_dir.as_ref().to_path_buf()),
            tables: HashMap::new(),
        }
    }

    /// Create an empty manager (no table directory)
    pub fn empty() -> Self {
        Bc5dTableManager {
            table_dir: None,
            tables: HashMap::new(),
        }
    }

    /// Get or load the table for a caliber
    ///
    /// Tables are cached after first load.
    pub fn get_table(&mut self, caliber: f64) -> Result<&Bc5dTable, Bc5dError> {
        let caliber_key = caliber_to_key(caliber);

        // Check if already loaded
        if self.tables.contains_key(&caliber_key) {
            return Ok(self.tables.get(&caliber_key).unwrap());
        }

        // Need to load
        let table_dir = self.table_dir.as_ref().ok_or(Bc5dError::NoTableDirectory)?;
        let table_path = find_table_file(table_dir, caliber)?;
        let table = Bc5dTable::load(&table_path)?;
        self.tables.insert(caliber_key, table);
        Ok(self.tables.get(&caliber_key).unwrap())
    }

    /// Look up BC correction for a bullet
    pub fn lookup(
        &mut self,
        caliber: f64,
        weight_grains: f64,
        base_bc: f64,
        muzzle_velocity: f64,
        current_velocity: f64,
        drag_type: &str,
    ) -> Result<f64, Bc5dError> {
        let table = self.get_table(caliber)?;
        Ok(table.lookup(weight_grains, base_bc, muzzle_velocity, current_velocity, drag_type))
    }

    /// Get effective BC with correction applied
    pub fn get_effective_bc(
        &mut self,
        caliber: f64,
        weight_grains: f64,
        base_bc: f64,
        muzzle_velocity: f64,
        current_velocity: f64,
        drag_type: &str,
    ) -> Result<f64, Bc5dError> {
        let table = self.get_table(caliber)?;
        Ok(table.get_effective_bc(weight_grains, base_bc, muzzle_velocity, current_velocity, drag_type))
    }

    /// Check if a table is available for a caliber
    pub fn has_table(&self, caliber: f64) -> bool {
        if let Some(ref table_dir) = self.table_dir {
            find_table_file(table_dir, caliber).is_ok()
        } else {
            false
        }
    }

    /// List available calibers in the table directory
    pub fn available_calibers(&self) -> Vec<f64> {
        let mut calibers = Vec::new();
        if let Some(ref table_dir) = self.table_dir {
            if let Ok(entries) = std::fs::read_dir(table_dir) {
                for entry in entries.flatten() {
                    let path = entry.path();
                    if let Some(ext) = path.extension() {
                        if ext == "bin" {
                            if let Some(stem) = path.file_stem() {
                                let name = stem.to_string_lossy();
                                if name.starts_with("bc5d_") {
                                    // Parse caliber from filename (e.g., bc5d_308.bin -> 0.308)
                                    if let Ok(cal_int) = name[5..].parse::<i32>() {
                                        calibers.push(cal_int as f64 / 1000.0);
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        calibers.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
        calibers
    }
}

/// Convert caliber to integer key (multiply by 1000)
fn caliber_to_key(caliber: f64) -> i32 {
    (caliber * 1000.0).round() as i32
}

/// Find the table file for a caliber
fn find_table_file(table_dir: &Path, caliber: f64) -> Result<PathBuf, Bc5dError> {
    let caliber_int = (caliber * 1000.0).round() as i32;
    let filename = format!("bc5d_{}.bin", caliber_int);
    let path = table_dir.join(&filename);

    if path.exists() {
        return Ok(path);
    }

    // Try common variations
    let variations = [
        format!("bc5d_{:03}.bin", caliber_int),
        format!("bc5d_0{}.bin", caliber_int),
    ];

    for var in &variations {
        let var_path = table_dir.join(var);
        if var_path.exists() {
            return Ok(var_path);
        }
    }

    Err(Bc5dError::TableNotFound(caliber))
}

// Helper functions for reading binary data

fn read_u32<R: Read>(reader: &mut R) -> Result<u32, std::io::Error> {
    let mut buf = [0u8; 4];
    reader.read_exact(&mut buf)?;
    Ok(u32::from_le_bytes(buf))
}

fn read_u64<R: Read>(reader: &mut R) -> Result<u64, std::io::Error> {
    let mut buf = [0u8; 8];
    reader.read_exact(&mut buf)?;
    Ok(u64::from_le_bytes(buf))
}

fn read_f32<R: Read>(reader: &mut R) -> Result<f32, std::io::Error> {
    let mut buf = [0u8; 4];
    reader.read_exact(&mut buf)?;
    Ok(f32::from_le_bytes(buf))
}

fn read_f32_array<R: Read>(reader: &mut R, count: usize) -> Result<Vec<f32>, std::io::Error> {
    // Defensive bounds: reject absurd lengths from corrupt/hostile files before
    // allocating, and guard the byte-count multiply against overflow.
    const MAX_ELEMS: usize = 64_000_000; // 256 MB of f32
    if count > MAX_ELEMS {
        return Err(std::io::Error::new(
            std::io::ErrorKind::InvalidData,
            "f32 array length too large",
        ));
    }
    let byte_len = count.checked_mul(4).ok_or_else(|| {
        std::io::Error::new(std::io::ErrorKind::InvalidData, "f32 array length overflow")
    })?;
    let mut data = vec![0f32; count];
    let mut buf = vec![0u8; byte_len];
    reader.read_exact(&mut buf)?;

    for (i, chunk) in buf.chunks_exact(4).enumerate() {
        data[i] = f32::from_le_bytes([chunk[0], chunk[1], chunk[2], chunk[3]]);
    }

    Ok(data)
}

/// Simple CRC32 (IEEE polynomial) implementation
pub(crate) fn crc32_ieee(data: &[u8]) -> u32 {
    const TABLE: [u32; 256] = make_crc32_table();
    let mut crc = 0xFFFFFFFFu32;
    for &byte in data {
        let idx = ((crc ^ byte as u32) & 0xFF) as usize;
        crc = (crc >> 8) ^ TABLE[idx];
    }
    !crc
}

const fn make_crc32_table() -> [u32; 256] {
    const POLY: u32 = 0xEDB88320;
    let mut table = [0u32; 256];
    let mut i = 0;
    while i < 256 {
        let mut crc = i as u32;
        let mut j = 0;
        while j < 8 {
            if crc & 1 != 0 {
                crc = (crc >> 1) ^ POLY;
            } else {
                crc >>= 1;
            }
            j += 1;
        }
        table[i] = crc;
        i += 1;
    }
    table
}

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

    fn create_test_table() -> Bc5dTable {
        // Create a small test table with known values
        let weight_bins = vec![100.0, 150.0, 200.0];
        let bc_bins = vec![0.3, 0.4, 0.5];
        let muzzle_vel_bins = vec![2500.0, 3000.0];
        let current_vel_bins = vec![1000.0, 2000.0, 3000.0];
        let num_drag_types = 2;

        // Total cells: 2 * 3 * 3 * 2 * 3 = 108
        let total = num_drag_types * weight_bins.len() * bc_bins.len() * muzzle_vel_bins.len() * current_vel_bins.len();
        let mut data = vec![1.0f32; total];

        // Set some non-uniform values for testing interpolation
        // At weight=150, bc=0.4, muzzle=2750 (interpolated), current=2000, G1
        // We'll set corners to test 4D interpolation
        data[0] = 0.95; // First corner
        data[total - 1] = 1.05; // Last corner

        Bc5dTable {
            caliber: 0.308,
            data,
            weight_bins,
            bc_bins,
            muzzle_vel_bins,
            current_vel_bins,
            num_drag_types,
            version: 2,
            api_version: "test".to_string(),
            timestamp: 0,
        }
    }

    fn create_single_cell_test_table() -> Bc5dTable {
        Bc5dTable {
            caliber: 0.308,
            data: vec![0.875],
            weight_bins: vec![168.0],
            bc_bins: vec![0.4],
            muzzle_vel_bins: vec![2500.0],
            current_vel_bins: vec![2000.0],
            num_drag_types: 1,
            version: 2,
            api_version: "test".to_string(),
            timestamp: 0,
        }
    }

    /// Serialize a table into the BC5D v2 `.bin` byte layout so we can exercise
    /// the `from_bytes` parser without depending on an external file.
    fn serialize_test_table(t: &Bc5dTable) -> Vec<u8> {
        let mut out = Vec::new();
        out.extend_from_slice(MAGIC);
        out.extend_from_slice(&t.version.to_le_bytes());
        out.extend_from_slice(&t.caliber.to_le_bytes());
        out.extend_from_slice(&0u32.to_le_bytes()); // flags
        out.extend_from_slice(&0u32.to_le_bytes()); // padding
        out.extend_from_slice(&(t.weight_bins.len() as u32).to_le_bytes());
        out.extend_from_slice(&(t.bc_bins.len() as u32).to_le_bytes());
        out.extend_from_slice(&(t.muzzle_vel_bins.len() as u32).to_le_bytes());
        out.extend_from_slice(&(t.current_vel_bins.len() as u32).to_le_bytes());
        out.extend_from_slice(&(t.num_drag_types as u32).to_le_bytes());
        out.extend_from_slice(&t.timestamp.to_le_bytes());

        // Checksum is CRC32 of bins + data, in declaration order.
        let mut checksum_data = Vec::new();
        for v in t.weight_bins.iter().chain(&t.bc_bins).chain(&t.muzzle_vel_bins)
            .chain(&t.current_vel_bins).chain(&t.data) {
            checksum_data.extend_from_slice(&v.to_le_bytes());
        }
        out.extend_from_slice(&crc32_ieee(&checksum_data).to_le_bytes());

        let mut api = [0u8; 16];
        let bytes = t.api_version.as_bytes();
        api[..bytes.len().min(16)].copy_from_slice(&bytes[..bytes.len().min(16)]);
        out.extend_from_slice(&api);
        out.extend_from_slice(&[0u8; 12]); // reserved

        for v in t.weight_bins.iter().chain(&t.bc_bins).chain(&t.muzzle_vel_bins)
            .chain(&t.current_vel_bins).chain(&t.data) {
            out.extend_from_slice(&v.to_le_bytes());
        }
        out
    }

    #[test]
    fn test_from_bytes_roundtrip() {
        let original = create_test_table();
        let bytes = serialize_test_table(&original);
        let parsed = Bc5dTable::from_bytes(&bytes).expect("from_bytes should parse");

        assert_eq!(parsed.caliber, original.caliber);
        assert_eq!(parsed.num_drag_types, original.num_drag_types);
        assert_eq!(parsed.weight_bins, original.weight_bins);
        assert_eq!(parsed.current_vel_bins, original.current_vel_bins);
        assert_eq!(parsed.data, original.data);
        assert_eq!(parsed.api_version, original.api_version);

        // A corrupted body must be rejected by the CRC check.
        let mut bad = bytes.clone();
        *bad.last_mut().unwrap() ^= 0xFF;
        assert!(Bc5dTable::from_bytes(&bad).is_err());
    }

    #[test]
    fn test_generate_segments() {
        // A table whose corrections are all exactly 1.0 carries no useful
        // correction, so generate_segments returns None (leave published BC).
        let mut uniform = create_test_table();
        uniform.data.iter_mut().for_each(|v| *v = 1.0);
        assert!(uniform
            .generate_segments(0.4, "G1", 150.0, Some(2700.0))
            .is_none());

        // A table with a real (0.9) correction across the sampled slice must
        // produce contiguous, descending velocity segments carrying bc*corr.
        let mut corrected = create_test_table();
        corrected.data.iter_mut().for_each(|v| *v = 0.9);
        let segments = corrected
            .generate_segments(0.4, "G1", 150.0, Some(2700.0))
            .expect("segments expected for a table with corrections");
        assert!(!segments.is_empty());
        for w in segments.windows(2) {
            // Bands are contiguous and descend in velocity.
            assert!((segments[0].velocity_max - w[0].velocity_max).abs() >= 0.0);
            assert!(w[0].velocity_min >= w[1].velocity_max - 1e-6);
        }
        for s in &segments {
            assert!((s.bc_value - 0.4 * 0.9).abs() < 1e-6); // base_bc * correction
            assert!(s.velocity_max > s.velocity_min);
        }
    }

    #[test]
    fn segment_schedule_carries_muzzle_corrected_fallback_bc() {
        let table = create_single_cell_test_table();
        let base_bc = 0.4;
        let schedule = table
            .generate_segment_schedule(base_bc, "G1", 168.0, 2500.0)
            .expect("uniform non-neutral correction should produce a schedule");
        let expected_fallback = table.get_effective_bc(168.0, base_bc, 2500.0, 2500.0, "G1");

        assert!(!schedule.segments.is_empty());
        assert_eq!(expected_fallback.to_bits(), (base_bc * 0.875).to_bits());
        assert_eq!(schedule.fallback_bc.to_bits(), expected_fallback.to_bits());
    }

    #[test]
    fn test_interp_idx_in_range() {
        let table = create_test_table();

        // Test middle of range
        let (idx, weight) = table.interp_idx(125.0, &table.weight_bins);
        assert_eq!(idx, 0);
        assert!((weight - 0.5).abs() < 0.01);

        // Test at bin boundary
        let (idx, weight) = table.interp_idx(150.0, &table.weight_bins);
        assert_eq!(idx, 0);
        assert!((weight - 1.0).abs() < 0.01);
    }

    #[test]
    fn test_interp_idx_out_of_range() {
        let table = create_test_table();

        // Test below range
        let (idx, weight) = table.interp_idx(50.0, &table.weight_bins);
        assert_eq!(idx, 0);
        assert_eq!(weight, 0.0);

        // Test above range
        let (idx, weight) = table.interp_idx(250.0, &table.weight_bins);
        assert_eq!(idx, 1); // len - 2
        assert_eq!(weight, 1.0);
    }

    #[test]
    fn test_interp_idx_nan_defaults_to_first_bin() {
        let table = create_test_table();

        assert_eq!(table.interp_idx(f32::NAN, &table.weight_bins), (0, 0.0));
        assert_eq!(
            table.interp_idx(f32::NEG_INFINITY, &table.weight_bins),
            (0, 0.0)
        );
        assert_eq!(
            table.interp_idx(f32::INFINITY, &table.weight_bins),
            (1, 1.0)
        );
    }

    #[test]
    fn test_lookup_nan_with_single_bin_axes_uses_only_cell() {
        let table = create_single_cell_test_table();

        assert_eq!(table.lookup(f64::NAN, 0.4, 2500.0, 2000.0, "G1"), 0.875);
    }

    #[test]
    fn test_lookup_non_finite_table_cells_are_neutral() {
        for cell in [f32::NAN, f32::INFINITY, f32::NEG_INFINITY] {
            let mut source = create_test_table();
            source.data.fill(cell);
            let bytes = serialize_test_table(&source);
            let table = Bc5dTable::from_bytes(&bytes).expect("CRC-valid table should load");

            assert_eq!(
                table.lookup(125.0, 0.35, 2750.0, 1500.0, "G1"),
                1.0,
                "non-finite cell {cell:?} must produce a neutral correction"
            );
            assert_eq!(
                table.get_effective_bc(125.0, 0.35, 2750.0, 1500.0, "G1"),
                0.35,
                "neutral correction must preserve base BC for {cell:?}"
            );
        }
    }

    #[test]
    fn test_lookup_returns_valid_range() {
        let table = create_test_table();

        let correction = table.lookup(150.0, 0.4, 2750.0, 2000.0, "G1");
        assert!(correction >= 0.5 && correction <= 1.5);

        let correction = table.lookup(150.0, 0.4, 2750.0, 2000.0, "G7");
        assert!(correction >= 0.5 && correction <= 1.5);
    }

    #[test]
    fn test_effective_bc() {
        let table = create_test_table();

        let base_bc = 0.4;
        let effective = table.get_effective_bc(150.0, base_bc, 2750.0, 2000.0, "G1");

        // Effective BC should be base_bc * correction
        assert!(effective >= base_bc * 0.5 && effective <= base_bc * 1.5);
    }

    #[test]
    fn test_caliber_to_key() {
        assert_eq!(caliber_to_key(0.308), 308);
        assert_eq!(caliber_to_key(0.224), 224);
        assert_eq!(caliber_to_key(0.338), 338);
    }

    #[test]
    fn test_table_metadata() {
        let table = create_test_table();
        assert!((table.caliber() - 0.308).abs() < 0.001);
        assert_eq!(table.version(), 2);
        assert_eq!(table.api_version(), "test");
    }

    #[test]
    fn test_crc32() {
        // Test with known CRC32 value
        let data = b"123456789";
        let crc = crc32_ieee(data);
        assert_eq!(crc, 0xCBF43926);
    }
}