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        // Estimate bullet length from caliber (typical rifle bullets are 3-4 calibers long)
244        let estimated_length = caliber_inches * 3.5;
245        self.lookup(bc, bc_type, mass_grains, estimated_length, velocity_fps)
246    }
247
248    /// Find interpolation index and weight for a value in bins
249    fn interp_idx(&self, value: f32, bins: &[f32]) -> (usize, f64) {
250        if bins.is_empty() {
251            return (0, 0.0);
252        }
253
254        // Handle out of range
255        if value <= bins[0] {
256            return (0, 0.0);
257        }
258        if value >= bins[bins.len() - 1] {
259            return (bins.len().saturating_sub(2), 1.0);
260        }
261
262        // Binary search for interval
263        let idx = match bins
264            .binary_search_by(|probe| probe.partial_cmp(&value).unwrap_or(std::cmp::Ordering::Equal))
265        {
266            Ok(i) => i.saturating_sub(1).min(bins.len() - 2),
267            Err(i) => i.saturating_sub(1).min(bins.len() - 2),
268        };
269
270        // Calculate interpolation weight
271        let low = bins[idx];
272        let high = bins[idx + 1];
273        let weight = if high > low {
274            ((value - low) / (high - low)) as f64
275        } else {
276            0.0
277        };
278
279        (idx, weight)
280    }
281
282    /// Calculate flat array index from 5D indices
283    fn flat_index(&self, type_idx: usize, bc_idx: usize, mass_idx: usize, length_idx: usize, vel_idx: usize) -> usize {
284        let n_bc = self.bc_bins.len();
285        let n_mass = self.mass_bins.len();
286        let n_length = self.length_bins.len();
287        let n_velocity = self.velocity_bins.len();
288
289        type_idx * (n_bc * n_mass * n_length * n_velocity)
290            + bc_idx * (n_mass * n_length * n_velocity)
291            + mass_idx * (n_length * n_velocity)
292            + length_idx * n_velocity
293            + vel_idx
294    }
295
296    /// Get table version
297    pub fn version(&self) -> u32 {
298        self.version
299    }
300
301    /// Get total number of cells in the table
302    pub fn total_cells(&self) -> usize {
303        self.data.len()
304    }
305
306    /// Get table dimensions as a string
307    pub fn dimensions_str(&self) -> String {
308        format!(
309            "{}x{}x{}x{}x{}",
310            self.bc_bins.len(),
311            self.mass_bins.len(),
312            self.length_bins.len(),
313            self.velocity_bins.len(),
314            self.num_types
315        )
316    }
317}
318
319// Helper functions for reading binary data
320
321fn read_u32<R: Read>(reader: &mut R) -> Result<u32, std::io::Error> {
322    let mut buf = [0u8; 4];
323    reader.read_exact(&mut buf)?;
324    Ok(u32::from_le_bytes(buf))
325}
326
327fn read_u64<R: Read>(reader: &mut R) -> Result<u64, std::io::Error> {
328    let mut buf = [0u8; 8];
329    reader.read_exact(&mut buf)?;
330    Ok(u64::from_le_bytes(buf))
331}
332
333fn read_f32_array<R: Read>(reader: &mut R, count: usize) -> Result<Vec<f32>, std::io::Error> {
334    // Defensive bounds: reject absurd lengths from corrupt/hostile files before allocating,
335    // and guard the byte-count multiply against overflow (mirrors bc_table_5d.rs).
336    const MAX_ELEMS: usize = 64_000_000; // 256 MB of f32
337    if count > MAX_ELEMS {
338        return Err(std::io::Error::new(
339            std::io::ErrorKind::InvalidData,
340            "f32 array length too large",
341        ));
342    }
343    let byte_len = count.checked_mul(4).ok_or_else(|| {
344        std::io::Error::new(std::io::ErrorKind::InvalidData, "f32 array length overflow")
345    })?;
346    let mut data = vec![0f32; count];
347    let mut buf = vec![0u8; byte_len];
348    reader.read_exact(&mut buf)?;
349
350    for (i, chunk) in buf.chunks_exact(4).enumerate() {
351        data[i] = f32::from_le_bytes([chunk[0], chunk[1], chunk[2], chunk[3]]);
352    }
353
354    Ok(data)
355}
356
357#[cfg(test)]
358mod tests {
359    use super::*;
360
361    #[test]
362    fn test_interp_idx_in_range() {
363        let table = BcCorrectionTable {
364            data: vec![1.0; 100],
365            bc_bins: vec![0.1, 0.2, 0.3, 0.4, 0.5],
366            mass_bins: vec![100.0, 150.0, 200.0],
367            length_bins: vec![1.0, 1.2, 1.4],
368            velocity_bins: vec![3000.0, 2500.0, 2000.0],
369            num_types: 2,
370            version: 1,
371        };
372
373        // Test middle of range
374        let (idx, weight) = table.interp_idx(0.25, &table.bc_bins);
375        assert_eq!(idx, 1);
376        assert!((weight - 0.5).abs() < 0.01);
377
378        // Test at bin boundary
379        let (idx, weight) = table.interp_idx(0.2, &table.bc_bins);
380        assert_eq!(idx, 0);
381        assert!((weight - 1.0).abs() < 0.01);
382    }
383
384    #[test]
385    fn test_interp_idx_out_of_range() {
386        let table = BcCorrectionTable {
387            data: vec![1.0; 100],
388            bc_bins: vec![0.1, 0.2, 0.3, 0.4, 0.5],
389            mass_bins: vec![100.0, 150.0, 200.0],
390            length_bins: vec![1.0, 1.2, 1.4],
391            velocity_bins: vec![3000.0, 2500.0, 2000.0],
392            num_types: 2,
393            version: 1,
394        };
395
396        // Test below range
397        let (idx, weight) = table.interp_idx(0.05, &table.bc_bins);
398        assert_eq!(idx, 0);
399        assert_eq!(weight, 0.0);
400
401        // Test above range
402        let (idx, weight) = table.interp_idx(0.6, &table.bc_bins);
403        assert_eq!(idx, 3); // len - 2
404        assert_eq!(weight, 1.0);
405    }
406}
ballistics_engine/bc_table.rs

ballistics_engine/
bc_table.rs
