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 _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        Ok(BcCorrectionTable {
143            data,
144            bc_bins,
145            mass_bins,
146            length_bins,
147            velocity_bins,
148            num_types,
149            version,
150        })
151    }
152
153    /// Look up a BC correction factor with 5D linear interpolation
154    ///
155    /// # Arguments
156    /// * `bc` - Published BC value
157    /// * `bc_type` - "G1" or "G7"
158    /// * `mass` - Bullet mass in grains
159    /// * `length` - Bullet length in inches
160    /// * `velocity` - Current velocity in fps
161    ///
162    /// # Returns
163    /// Correction factor (multiply published BC by this value)
164    pub fn lookup(&self, bc: f64, bc_type: &str, mass: f64, length: f64, velocity: f64) -> f64 {
165        // Get type index (0 = G1, 1 = G7), clamped to the table's type axis so a
166        // G1-only table (num_types == 1) queried with "G7" stays in bounds instead of
167        // producing an out-of-range flat_index (mirrors Bc5dTable::lookup).
168        let type_idx =
169            (if bc_type.to_uppercase() == "G1" { 0 } else { 1 }).min(self.num_types.saturating_sub(1));
170
171        // Find interpolation indices and weights for each dimension
172        let (bc_idx, bc_weight) = self.interp_idx(bc as f32, &self.bc_bins);
173        let (mass_idx, mass_weight) = self.interp_idx(mass as f32, &self.mass_bins);
174        let (length_idx, length_weight) = self.interp_idx(length as f32, &self.length_bins);
175        let (vel_idx, vel_weight) = self.interp_idx(velocity as f32, &self.velocity_bins);
176
177        // 4D linear interpolation (type is discrete, not interpolated)
178        let mut result = 0.0f64;
179
180        for di in 0..2 {
181            for dj in 0..2 {
182                for dk in 0..2 {
183                    for dl in 0..2 {
184                        // Calculate weight for this corner
185                        let w = (if di == 0 { 1.0 - bc_weight } else { bc_weight })
186                            * (if dj == 0 { 1.0 - mass_weight } else { mass_weight })
187                            * (if dk == 0 { 1.0 - length_weight } else { length_weight })
188                            * (if dl == 0 { 1.0 - vel_weight } else { vel_weight });
189
190                        // Get clamped indices
191                        let i = (bc_idx + di).min(self.bc_bins.len() - 1);
192                        let j = (mass_idx + dj).min(self.mass_bins.len() - 1);
193                        let k = (length_idx + dk).min(self.length_bins.len() - 1);
194                        let l = (vel_idx + dl).min(self.velocity_bins.len() - 1);
195
196                        // Calculate flat index
197                        let idx = self.flat_index(type_idx, i, j, k, l);
198                        // Bounds-safe access (true defense-in-depth): the continuous axes are
199                        // clamped above and type_idx is clamped to num_types, so idx is always
200                        // in bounds for any valid table; the unwrap_or only guards corrupt data.
201                        result += w * self.data.get(idx).copied().unwrap_or(0.0) as f64;
202                    }
203                }
204            }
205        }
206
207        result
208    }
209
210    /// Get the BC correction factor for a given bullet and velocity
211    /// This is a convenience method that estimates bullet length from mass and caliber
212    ///
213    /// # Arguments
214    /// * `bc` - Published BC value
215    /// * `bc_type` - "G1" or "G7"
216    /// * `mass_grains` - Bullet mass in grains
217    /// * `caliber_inches` - Bullet caliber in inches
218    /// * `velocity_fps` - Current velocity in fps
219    pub fn lookup_with_caliber(
220        &self,
221        bc: f64,
222        bc_type: &str,
223        mass_grains: f64,
224        caliber_inches: f64,
225        velocity_fps: f64,
226    ) -> f64 {
227        // Estimate bullet length from caliber (typical rifle bullets are 3-4 calibers long)
228        let estimated_length = caliber_inches * 3.5;
229        self.lookup(bc, bc_type, mass_grains, estimated_length, velocity_fps)
230    }
231
232    /// Find interpolation index and weight for a value in bins
233    fn interp_idx(&self, value: f32, bins: &[f32]) -> (usize, f64) {
234        if bins.is_empty() {
235            return (0, 0.0);
236        }
237
238        // Handle out of range
239        if value <= bins[0] {
240            return (0, 0.0);
241        }
242        if value >= bins[bins.len() - 1] {
243            return (bins.len().saturating_sub(2), 1.0);
244        }
245
246        // Binary search for interval
247        let idx = match bins
248            .binary_search_by(|probe| probe.partial_cmp(&value).unwrap_or(std::cmp::Ordering::Equal))
249        {
250            Ok(i) => i.saturating_sub(1).min(bins.len() - 2),
251            Err(i) => i.saturating_sub(1).min(bins.len() - 2),
252        };
253
254        // Calculate interpolation weight
255        let low = bins[idx];
256        let high = bins[idx + 1];
257        let weight = if high > low {
258            ((value - low) / (high - low)) as f64
259        } else {
260            0.0
261        };
262
263        (idx, weight)
264    }
265
266    /// Calculate flat array index from 5D indices
267    fn flat_index(&self, type_idx: usize, bc_idx: usize, mass_idx: usize, length_idx: usize, vel_idx: usize) -> usize {
268        let n_bc = self.bc_bins.len();
269        let n_mass = self.mass_bins.len();
270        let n_length = self.length_bins.len();
271        let n_velocity = self.velocity_bins.len();
272
273        type_idx * (n_bc * n_mass * n_length * n_velocity)
274            + bc_idx * (n_mass * n_length * n_velocity)
275            + mass_idx * (n_length * n_velocity)
276            + length_idx * n_velocity
277            + vel_idx
278    }
279
280    /// Get table version
281    pub fn version(&self) -> u32 {
282        self.version
283    }
284
285    /// Get total number of cells in the table
286    pub fn total_cells(&self) -> usize {
287        self.data.len()
288    }
289
290    /// Get table dimensions as a string
291    pub fn dimensions_str(&self) -> String {
292        format!(
293            "{}x{}x{}x{}x{}",
294            self.bc_bins.len(),
295            self.mass_bins.len(),
296            self.length_bins.len(),
297            self.velocity_bins.len(),
298            self.num_types
299        )
300    }
301}
302
303// Helper functions for reading binary data
304
305fn read_u32<R: Read>(reader: &mut R) -> Result<u32, std::io::Error> {
306    let mut buf = [0u8; 4];
307    reader.read_exact(&mut buf)?;
308    Ok(u32::from_le_bytes(buf))
309}
310
311fn read_u64<R: Read>(reader: &mut R) -> Result<u64, std::io::Error> {
312    let mut buf = [0u8; 8];
313    reader.read_exact(&mut buf)?;
314    Ok(u64::from_le_bytes(buf))
315}
316
317fn read_f32_array<R: Read>(reader: &mut R, count: usize) -> Result<Vec<f32>, std::io::Error> {
318    // Defensive bounds: reject absurd lengths from corrupt/hostile files before allocating,
319    // and guard the byte-count multiply against overflow (mirrors bc_table_5d.rs).
320    const MAX_ELEMS: usize = 64_000_000; // 256 MB of f32
321    if count > MAX_ELEMS {
322        return Err(std::io::Error::new(
323            std::io::ErrorKind::InvalidData,
324            "f32 array length too large",
325        ));
326    }
327    let byte_len = count.checked_mul(4).ok_or_else(|| {
328        std::io::Error::new(std::io::ErrorKind::InvalidData, "f32 array length overflow")
329    })?;
330    let mut data = vec![0f32; count];
331    let mut buf = vec![0u8; byte_len];
332    reader.read_exact(&mut buf)?;
333
334    for (i, chunk) in buf.chunks_exact(4).enumerate() {
335        data[i] = f32::from_le_bytes([chunk[0], chunk[1], chunk[2], chunk[3]]);
336    }
337
338    Ok(data)
339}
340
341#[cfg(test)]
342mod tests {
343    use super::*;
344
345    #[test]
346    fn test_interp_idx_in_range() {
347        let table = BcCorrectionTable {
348            data: vec![1.0; 100],
349            bc_bins: vec![0.1, 0.2, 0.3, 0.4, 0.5],
350            mass_bins: vec![100.0, 150.0, 200.0],
351            length_bins: vec![1.0, 1.2, 1.4],
352            velocity_bins: vec![3000.0, 2500.0, 2000.0],
353            num_types: 2,
354            version: 1,
355        };
356
357        // Test middle of range
358        let (idx, weight) = table.interp_idx(0.25, &table.bc_bins);
359        assert_eq!(idx, 1);
360        assert!((weight - 0.5).abs() < 0.01);
361
362        // Test at bin boundary
363        let (idx, weight) = table.interp_idx(0.2, &table.bc_bins);
364        assert_eq!(idx, 0);
365        assert!((weight - 1.0).abs() < 0.01);
366    }
367
368    #[test]
369    fn test_interp_idx_out_of_range() {
370        let table = BcCorrectionTable {
371            data: vec![1.0; 100],
372            bc_bins: vec![0.1, 0.2, 0.3, 0.4, 0.5],
373            mass_bins: vec![100.0, 150.0, 200.0],
374            length_bins: vec![1.0, 1.2, 1.4],
375            velocity_bins: vec![3000.0, 2500.0, 2000.0],
376            num_types: 2,
377            version: 1,
378        };
379
380        // Test below range
381        let (idx, weight) = table.interp_idx(0.05, &table.bc_bins);
382        assert_eq!(idx, 0);
383        assert_eq!(weight, 0.0);
384
385        // Test above range
386        let (idx, weight) = table.interp_idx(0.6, &table.bc_bins);
387        assert_eq!(idx, 3); // len - 2
388        assert_eq!(weight, 1.0);
389    }
390}
ballistics_engine/bc_table.rs

ballistics_engine/
bc_table.rs
