ballistics_engine/
bc_table.rs

1// BC Correction Table - 5D Lookup Table with Linear Interpolation
2//
3// This module provides offline BC corrections by loading a precomputed table
4// of correction factors. The table is indexed by:
5//   - BC value (published BC)
6//   - BC type (G1 or G7)
7//   - Bullet mass (grains)
8//   - Bullet length (inches)
9//   - Velocity (fps)
10//
11// Binary file format (BCCR):
12//   Header (64 bytes):
13//     - Magic: 4 bytes ('BCCR')
14//     - Version: 4 bytes (uint32)
15//     - Flags: 4 bytes (uint32)
16//     - num_bc: 4 bytes (uint32)
17//     - num_mass: 4 bytes (uint32)
18//     - num_length: 4 bytes (uint32)
19//     - num_velocity: 4 bytes (uint32)
20//     - num_bc_types: 4 bytes (uint32)
21//     - timestamp: 8 bytes (uint64)
22//     - checksum: 4 bytes (uint32, CRC32)
23//     - reserved: 16 bytes
24//   Bin definitions:
25//     - BC bins: num_bc * 4 bytes (float32)
26//     - Mass bins: num_mass * 4 bytes (float32)
27//     - Length bins: num_length * 4 bytes (float32)
28//     - Velocity bins: num_velocity * 4 bytes (float32)
29//   Data section:
30//     - Correction factors: total_cells * 4 bytes (float32)
31//     - Layout: [bc_type][bc][mass][length][velocity]
32
33use std::fs::File;
34use std::io::{BufReader, Read};
35use std::path::Path;
36
37/// BC correction table with 5D interpolation
38#[derive(Debug)]
39pub struct BcCorrectionTable {
40    /// Table data: [bc_type][bc][mass][length][velocity]
41    data: Vec<f32>,
42    /// BC bin values
43    bc_bins: Vec<f32>,
44    /// Mass bin values (grains)
45    mass_bins: Vec<f32>,
46    /// Length bin values (inches)
47    length_bins: Vec<f32>,
48    /// Velocity bin values (fps)
49    velocity_bins: Vec<f32>,
50    /// Number of BC types (typically 2: G1, G7)
51    num_types: usize,
52    /// Table version
53    version: u32,
54}
55
56/// Error type for BC table operations
57#[derive(Debug)]
58pub enum BcTableError {
59    IoError(std::io::Error),
60    InvalidMagic,
61    UnsupportedVersion(u32),
62    ChecksumMismatch,
63    InvalidDimensions,
64}
65
66impl std::fmt::Display for BcTableError {
67    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
68        match self {
69            BcTableError::IoError(e) => write!(f, "IO error: {}", e),
70            BcTableError::InvalidMagic => write!(f, "Invalid file magic (expected 'BCCR')"),
71            BcTableError::UnsupportedVersion(v) => write!(f, "Unsupported table version: {}", v),
72            BcTableError::ChecksumMismatch => write!(f, "Data checksum mismatch"),
73            BcTableError::InvalidDimensions => write!(f, "Invalid table dimensions"),
74        }
75    }
76}
77
78impl std::error::Error for BcTableError {}
79
80impl From<std::io::Error> for BcTableError {
81    fn from(e: std::io::Error) -> Self {
82        BcTableError::IoError(e)
83    }
84}
85
86impl BcCorrectionTable {
87    /// Load a BC correction table from a binary file
88    pub fn load<P: AsRef<Path>>(path: P) -> Result<Self, BcTableError> {
89        let file = File::open(path)?;
90        let mut reader = BufReader::new(file);
91
92        // Read header (64 bytes)
93        let mut magic = [0u8; 4];
94        reader.read_exact(&mut magic)?;
95        if &magic != b"BCCR" {
96            return Err(BcTableError::InvalidMagic);
97        }
98
99        let version = read_u32(&mut reader)?;
100        if version != 1 {
101            return Err(BcTableError::UnsupportedVersion(version));
102        }
103
104        let _flags = read_u32(&mut reader)?;
105        let num_bc = read_u32(&mut reader)? as usize;
106        let num_mass = read_u32(&mut reader)? as usize;
107        let num_length = read_u32(&mut reader)? as usize;
108        let num_velocity = read_u32(&mut reader)? as usize;
109        let num_types = read_u32(&mut reader)? as usize;
110        let _timestamp = read_u64(&mut reader)?;
111        let stored_checksum = read_u32(&mut reader)?;
112
113        // Skip reserved bytes
114        let mut reserved = [0u8; 16];
115        reader.read_exact(&mut reserved)?;
116
117        // Validate dimensions
118        if num_bc == 0 || num_mass == 0 || num_length == 0 || num_velocity == 0 || num_types == 0 {
119            return Err(BcTableError::InvalidDimensions);
120        }
121
122        // Read bin definitions
123        let bc_bins = read_f32_array(&mut reader, num_bc)?;
124        let mass_bins = read_f32_array(&mut reader, num_mass)?;
125        let length_bins = read_f32_array(&mut reader, num_length)?;
126        let velocity_bins = read_f32_array(&mut reader, num_velocity)?;
127
128        // Read data section. Bound the product with checked arithmetic so a corrupt/hostile
129        // file cannot overflow (debug panic / release wrap to a too-short `data`, which then
130        // panics on the unchecked index in lookup()) or trigger a huge OOM allocation.
131        // Mirrors the bc_table_5d.rs hardening.
132        const MAX_TOTAL_CELLS: usize = 64_000_000; // ~256 MB of f32; far above any real table
133        let total_cells = num_types
134            .checked_mul(num_bc)
135            .and_then(|x| x.checked_mul(num_mass))
136            .and_then(|x| x.checked_mul(num_length))
137            .and_then(|x| x.checked_mul(num_velocity))
138            .filter(|&n| n <= MAX_TOTAL_CELLS)
139            .ok_or(BcTableError::InvalidDimensions)?;
140        let data = read_f32_array(&mut reader, total_cells)?;
141
142        // MBA-953: verify the stored CRC32 (previously read and discarded, so corrupt-but-in-range
143        // files loaded silently). BCCR CRCs the DATA section ONLY — the f32 cells — unlike
144        // bc_table_5d which CRCs bins+data. Gate on a non-zero stored value so older files written
145        // without a checksum (stored == 0) still load. The f32 round-trip is bit-exact, so this
146        // reproduces the on-disk data bytes.
147        if stored_checksum != 0 {
148            let mut data_bytes = Vec::with_capacity(data.len() * 4);
149            for &v in &data {
150                data_bytes.extend_from_slice(&v.to_le_bytes());
151            }
152            let calculated = crate::bc_table_5d::crc32_ieee(&data_bytes);
153            if calculated != stored_checksum {
154                return Err(BcTableError::ChecksumMismatch);
155            }
156        }
157
158        Ok(BcCorrectionTable {
159            data,
160            bc_bins,
161            mass_bins,
162            length_bins,
163            velocity_bins,
164            num_types,
165            version,
166        })
167    }
168
169    /// Look up a BC correction factor with 5D linear interpolation
170    ///
171    /// # Arguments
172    /// * `bc` - Published BC value
173    /// * `bc_type` - "G1" or "G7"
174    /// * `mass` - Bullet mass in grains
175    /// * `length` - Bullet length in inches
176    /// * `velocity` - Current velocity in fps
177    ///
178    /// # Returns
179    /// Correction factor (multiply published BC by this value)
180    pub fn lookup(&self, bc: f64, bc_type: &str, mass: f64, length: f64, velocity: f64) -> f64 {
181        // Get type index (0 = G1, 1 = G7), clamped to the table's type axis so a
182        // G1-only table (num_types == 1) queried with "G7" stays in bounds instead of
183        // producing an out-of-range flat_index (mirrors Bc5dTable::lookup).
184        let type_idx =
185            (if bc_type.to_uppercase() == "G1" { 0 } else { 1 }).min(self.num_types.saturating_sub(1));
186
187        // Find interpolation indices and weights for each dimension
188        let (bc_idx, bc_weight) = self.interp_idx(bc as f32, &self.bc_bins);
189        let (mass_idx, mass_weight) = self.interp_idx(mass as f32, &self.mass_bins);
190        let (length_idx, length_weight) = self.interp_idx(length as f32, &self.length_bins);
191        let (vel_idx, vel_weight) = self.interp_idx(velocity as f32, &self.velocity_bins);
192
193        // 4D linear interpolation (type is discrete, not interpolated)
194        let mut result = 0.0f64;
195
196        for di in 0..2 {
197            for dj in 0..2 {
198                for dk in 0..2 {
199                    for dl in 0..2 {
200                        // Calculate weight for this corner
201                        let w = (if di == 0 { 1.0 - bc_weight } else { bc_weight })
202                            * (if dj == 0 { 1.0 - mass_weight } else { mass_weight })
203                            * (if dk == 0 { 1.0 - length_weight } else { length_weight })
204                            * (if dl == 0 { 1.0 - vel_weight } else { vel_weight });
205
206                        // Get clamped indices
207                        let i = (bc_idx + di).min(self.bc_bins.len() - 1);
208                        let j = (mass_idx + dj).min(self.mass_bins.len() - 1);
209                        let k = (length_idx + dk).min(self.length_bins.len() - 1);
210                        let l = (vel_idx + dl).min(self.velocity_bins.len() - 1);
211
212                        // Calculate flat index
213                        let idx = self.flat_index(type_idx, i, j, k, l);
214                        // Bounds-safe access (true defense-in-depth): the continuous axes are
215                        // clamped above and type_idx is clamped to num_types, so idx is always
216                        // in bounds for any valid table; the unwrap_or only guards corrupt data.
217                        result += w * self.data.get(idx).copied().unwrap_or(0.0) as f64;
218                    }
219                }
220            }
221        }
222
223        result
224    }
225
226    /// Get the BC correction factor for a given bullet and velocity
227    /// This is a convenience method that estimates bullet length from mass and caliber
228    ///
229    /// # Arguments
230    /// * `bc` - Published BC value
231    /// * `bc_type` - "G1" or "G7"
232    /// * `mass_grains` - Bullet mass in grains
233    /// * `caliber_inches` - Bullet caliber in inches
234    /// * `velocity_fps` - Current velocity in fps
235    pub fn lookup_with_caliber(
236        &self,
237        bc: f64,
238        bc_type: &str,
239        mass_grains: f64,
240        caliber_inches: f64,
241        velocity_fps: f64,
242    ) -> f64 {
243        // MBA-1135: estimate bullet length from mass + caliber (was a mass-blind caliber*3.5
244        // heuristic). Length is a real interpolation axis of this table, so a weight-aware
245        // estimate improves the correction. Fall back to the old heuristic if mass is unavailable.
246        let est_m = crate::stability::estimate_bullet_length_m(
247            caliber_inches * 0.0254,
248            mass_grains * 0.00006479891,
249        );
250        let estimated_length = if est_m > 0.0 {
251            est_m / 0.0254
252        } else {
253            caliber_inches * 3.5
254        };
255        self.lookup(bc, bc_type, mass_grains, estimated_length, velocity_fps)
256    }
257
258    /// Find interpolation index and weight for a value in bins
259    fn interp_idx(&self, value: f32, bins: &[f32]) -> (usize, f64) {
260        if bins.is_empty() {
261            return (0, 0.0);
262        }
263
264        // Handle out of range
265        if value <= bins[0] {
266            return (0, 0.0);
267        }
268        if value >= bins[bins.len() - 1] {
269            return (bins.len().saturating_sub(2), 1.0);
270        }
271
272        // Binary search for interval
273        let idx = match bins
274            .binary_search_by(|probe| probe.partial_cmp(&value).unwrap_or(std::cmp::Ordering::Equal))
275        {
276            Ok(i) => i.saturating_sub(1).min(bins.len() - 2),
277            Err(i) => i.saturating_sub(1).min(bins.len() - 2),
278        };
279
280        // Calculate interpolation weight
281        let low = bins[idx];
282        let high = bins[idx + 1];
283        let weight = if high > low {
284            ((value - low) / (high - low)) as f64
285        } else {
286            0.0
287        };
288
289        (idx, weight)
290    }
291
292    /// Calculate flat array index from 5D indices
293    fn flat_index(&self, type_idx: usize, bc_idx: usize, mass_idx: usize, length_idx: usize, vel_idx: usize) -> usize {
294        let n_bc = self.bc_bins.len();
295        let n_mass = self.mass_bins.len();
296        let n_length = self.length_bins.len();
297        let n_velocity = self.velocity_bins.len();
298
299        type_idx * (n_bc * n_mass * n_length * n_velocity)
300            + bc_idx * (n_mass * n_length * n_velocity)
301            + mass_idx * (n_length * n_velocity)
302            + length_idx * n_velocity
303            + vel_idx
304    }
305
306    /// Get table version
307    pub fn version(&self) -> u32 {
308        self.version
309    }
310
311    /// Get total number of cells in the table
312    pub fn total_cells(&self) -> usize {
313        self.data.len()
314    }
315
316    /// Get table dimensions as a string
317    pub fn dimensions_str(&self) -> String {
318        format!(
319            "{}x{}x{}x{}x{}",
320            self.bc_bins.len(),
321            self.mass_bins.len(),
322            self.length_bins.len(),
323            self.velocity_bins.len(),
324            self.num_types
325        )
326    }
327}
328
329// Helper functions for reading binary data
330
331fn read_u32<R: Read>(reader: &mut R) -> Result<u32, std::io::Error> {
332    let mut buf = [0u8; 4];
333    reader.read_exact(&mut buf)?;
334    Ok(u32::from_le_bytes(buf))
335}
336
337fn read_u64<R: Read>(reader: &mut R) -> Result<u64, std::io::Error> {
338    let mut buf = [0u8; 8];
339    reader.read_exact(&mut buf)?;
340    Ok(u64::from_le_bytes(buf))
341}
342
343fn read_f32_array<R: Read>(reader: &mut R, count: usize) -> Result<Vec<f32>, std::io::Error> {
344    // Defensive bounds: reject absurd lengths from corrupt/hostile files before allocating,
345    // and guard the byte-count multiply against overflow (mirrors bc_table_5d.rs).
346    const MAX_ELEMS: usize = 64_000_000; // 256 MB of f32
347    if count > MAX_ELEMS {
348        return Err(std::io::Error::new(
349            std::io::ErrorKind::InvalidData,
350            "f32 array length too large",
351        ));
352    }
353    let byte_len = count.checked_mul(4).ok_or_else(|| {
354        std::io::Error::new(std::io::ErrorKind::InvalidData, "f32 array length overflow")
355    })?;
356    let mut data = vec![0f32; count];
357    let mut buf = vec![0u8; byte_len];
358    reader.read_exact(&mut buf)?;
359
360    for (i, chunk) in buf.chunks_exact(4).enumerate() {
361        data[i] = f32::from_le_bytes([chunk[0], chunk[1], chunk[2], chunk[3]]);
362    }
363
364    Ok(data)
365}
366
367#[cfg(test)]
368mod tests {
369    use super::*;
370
371    #[test]
372    fn test_interp_idx_in_range() {
373        let table = BcCorrectionTable {
374            data: vec![1.0; 100],
375            bc_bins: vec![0.1, 0.2, 0.3, 0.4, 0.5],
376            mass_bins: vec![100.0, 150.0, 200.0],
377            length_bins: vec![1.0, 1.2, 1.4],
378            velocity_bins: vec![3000.0, 2500.0, 2000.0],
379            num_types: 2,
380            version: 1,
381        };
382
383        // Test middle of range
384        let (idx, weight) = table.interp_idx(0.25, &table.bc_bins);
385        assert_eq!(idx, 1);
386        assert!((weight - 0.5).abs() < 0.01);
387
388        // Test at bin boundary
389        let (idx, weight) = table.interp_idx(0.2, &table.bc_bins);
390        assert_eq!(idx, 0);
391        assert!((weight - 1.0).abs() < 0.01);
392    }
393
394    #[test]
395    fn test_interp_idx_out_of_range() {
396        let table = BcCorrectionTable {
397            data: vec![1.0; 100],
398            bc_bins: vec![0.1, 0.2, 0.3, 0.4, 0.5],
399            mass_bins: vec![100.0, 150.0, 200.0],
400            length_bins: vec![1.0, 1.2, 1.4],
401            velocity_bins: vec![3000.0, 2500.0, 2000.0],
402            num_types: 2,
403            version: 1,
404        };
405
406        // Test below range
407        let (idx, weight) = table.interp_idx(0.05, &table.bc_bins);
408        assert_eq!(idx, 0);
409        assert_eq!(weight, 0.0);
410
411        // Test above range
412        let (idx, weight) = table.interp_idx(0.6, &table.bc_bins);
413        assert_eq!(idx, 3); // len - 2
414        assert_eq!(weight, 1.0);
415    }
416}
ballistics_engine/bc_table.rs

ballistics_engine/
bc_table.rs
