Skip to main content

ballistics_engine/
wind.rs

1use nalgebra::Vector3;
2use std::cmp::Ordering;
3use std::f64::consts::PI;
4
5/// Conversion constant from KMH to MPS
6const KMH_TO_MPS: f64 = 1000.0 / 3600.0;
7
8/// Wind segment: (speed_kmh, angle_deg, until_distance_m)
9/// This matches the Python WindSock interface
10pub type WindSegment = (f64, f64, f64);
11
12/// Sort wind segments by their `until_distance_m` threshold.
13///
14/// Shared by [`WindSock`] and the low-level trajectory integrator so every segmented-wind path
15/// applies the same interval ordering.
16pub(crate) fn sort_wind_segments_by_distance(segments: &mut [WindSegment]) {
17    segments.sort_by(|a, b| match (a.2.is_nan(), b.2.is_nan()) {
18        (true, true) => Ordering::Equal,
19        (true, false) => Ordering::Greater,
20        (false, true) => Ordering::Less,
21        (false, false) => {
22            a.2.partial_cmp(&b.2)
23                .expect("non-NaN distances are ordered")
24        }
25    });
26}
27
28/// Wind condition handler for trajectory calculations
29#[derive(Debug, Clone)]
30pub struct WindSock {
31    /// Sorted wind segments by distance
32    winds: Vec<WindSegment>,
33    /// Precomputed wind vector for each segment (parallel to `winds`). The Monte-Carlo RK4
34    /// kernel queries wind 4x per step, so caching avoids recomputing sin/cos every call.
35    wind_vecs: Vec<Vector3<f64>>,
36    /// Current segment index
37    current: usize,
38    /// Distance where next segment starts
39    next_range: f64,
40    /// Current wind vector
41    current_vec: Vector3<f64>,
42}
43
44impl WindSock {
45    /// Create a new WindSock from wind segments
46    ///
47    /// Args:
48    ///     segments: List of (speed_kmh, angle_deg, until_distance_m) tuples
49    pub fn new(mut segments: Vec<WindSegment>) -> Self {
50        // Sort segments by distance, handling NaN safely by treating it as greater than any value
51        sort_wind_segments_by_distance(&mut segments);
52
53        // Precompute each segment's wind vector once (depends only on its speed/angle).
54        let wind_vecs: Vec<Vector3<f64>> = segments.iter().map(Self::calc_vec).collect();
55
56        let (current, next_range, current_vec) = if segments.is_empty() {
57            (0, f64::INFINITY, Vector3::zeros())
58        } else {
59            (0, segments[0].2, wind_vecs[0])
60        };
61
62        WindSock {
63            winds: segments,
64            wind_vecs,
65            current,
66            next_range,
67            current_vec,
68        }
69    }
70
71    /// Calculate wind vector from wind segment
72    fn calc_vec(seg: &WindSegment) -> Vector3<f64> {
73        let (speed_kmh, angle_deg, _) = *seg;
74
75        // Convert kmh to m/s
76        let speed_mps = speed_kmh * KMH_TO_MPS;
77        let angle_rad = angle_deg * PI / 180.0;
78
79        // Wind convention (matching trajectory coordinates):
80        // 0° = headwind (from front, affects -x downrange)
81        // 90° = wind from right (affects -z lateral)
82        // 180° = tailwind (from back, affects +x downrange)
83        // 270° = wind from left (affects +z lateral)
84        //
85        // McCoy convention: x=downrange, y=vertical, z=lateral
86        Vector3::new(
87            -speed_mps * angle_rad.cos(), // x (downrange - head/tail component)
88            0.0,                          // y (vertical)
89            -speed_mps * angle_rad.sin(), // z (lateral - crosswind component)
90        )
91    }
92
93    /// Get wind vector for a given range
94    ///
95    /// Note: This modifies internal state and expects monotonically increasing ranges
96    /// For trajectory integration, we need a stateless version
97    pub fn vector_for_range(&mut self, range_m: f64) -> Vector3<f64> {
98        // Handle NaN
99        if range_m.is_nan() {
100            return Vector3::zeros();
101        }
102
103        // Advance the cursor across however many segments the query skipped (a single `if`
104        // returned a stale vector when a monotonic query jumped past a whole short segment).
105        while range_m >= self.next_range && self.current < self.winds.len() {
106            self.current += 1;
107            if self.current >= self.winds.len() {
108                self.current_vec = Vector3::zeros();
109                self.next_range = f64::INFINITY;
110            } else {
111                self.current_vec = self.wind_vecs[self.current];
112                self.next_range = self.winds[self.current].2;
113            }
114        }
115
116        self.current_vec
117    }
118
119    /// Get wind vector for a given range (stateless version)
120    ///
121    /// This version doesn't modify internal state and is safe for numerical integration
122    /// where the same range might be queried multiple times or out of order
123    pub fn vector_for_range_stateless(&self, range_m: f64) -> Vector3<f64> {
124        // Handle NaN
125        if range_m.is_nan() {
126            return Vector3::zeros();
127        }
128
129        // Find the appropriate segment (precomputed vector — no per-call trig).
130        for (i, segment) in self.winds.iter().enumerate() {
131            if range_m < segment.2 {
132                return self.wind_vecs[i];
133            }
134        }
135
136        // Beyond all segments
137        Vector3::zeros()
138    }
139}
140
141/// Parse a `"SPEED:ANGLE:UNTIL_DISTANCE"` string into a [`WindSegment`]
142/// `(speed_kmh, angle_deg, until_distance_m)`.
143///
144/// `imperial`: when true, SPEED is mph and UNTIL_DISTANCE is yards; otherwise
145/// SPEED is m/s and UNTIL_DISTANCE is meters. ANGLE is always degrees in the
146/// wind-FROM convention (0 = headwind, 90 = from the right). Shared by the CLI
147/// (`--wind-segment`) and the WASM front-ends so they parse identically.
148pub fn parse_wind_segment_str(s: &str, imperial: bool) -> Result<WindSegment, String> {
149    let parts: Vec<&str> = s.split(':').collect();
150    if parts.len() != 3 {
151        return Err(format!(
152            "invalid wind segment '{s}': expected SPEED:ANGLE:UNTIL_DISTANCE (three colon-separated numbers)"
153        ));
154    }
155    let num = |i: usize, name: &str| -> Result<f64, String> {
156        parts[i].trim().parse::<f64>().map_err(|_| {
157            format!("invalid wind segment '{s}': {name} '{}' is not a number", parts[i])
158        })
159    };
160    let speed = num(0, "speed")?;
161    let angle = num(1, "angle")?;
162    let until = num(2, "until-distance")?;
163    if !speed.is_finite() || !angle.is_finite() || !until.is_finite() {
164        return Err(format!(
165            "invalid wind segment '{s}': speed, angle, and until-distance must be finite numbers"
166        ));
167    }
168    if speed < 0.0 {
169        return Err(format!("invalid wind segment '{s}': speed must be >= 0"));
170    }
171    if until <= 0.0 {
172        return Err(format!("invalid wind segment '{s}': until-distance must be > 0"));
173    }
174    let (speed_kmh, until_m) = if imperial {
175        (speed * 1.609344, until * 0.9144) // mph -> km/h, yards -> meters
176    } else {
177        (speed * 3.6, until) // m/s -> km/h, meters -> meters
178    };
179    Ok((speed_kmh, angle, until_m))
180}
181
182#[cfg(test)]
183mod tests {
184    use super::*;
185
186    #[test]
187    fn segment_sort_is_stable_and_places_nan_endpoints_last() {
188        let mut segments = vec![
189            (10.0, 0.0, f64::NAN),
190            (20.0, 0.0, 100.0),
191            (30.0, 0.0, 100.0),
192            (40.0, 0.0, f64::INFINITY),
193            (50.0, 0.0, f64::NEG_INFINITY),
194            (60.0, 0.0, f64::NAN),
195        ];
196
197        sort_wind_segments_by_distance(&mut segments);
198
199        assert_eq!(segments[0].0, 50.0); // -inf first
200        assert_eq!(segments[1].0, 20.0); // equal endpoints retain input order
201        assert_eq!(segments[2].0, 30.0);
202        assert_eq!(segments[3].0, 40.0); // +inf after finite endpoints
203        assert_eq!(segments[4].0, 10.0); // NaNs last and stable
204        assert_eq!(segments[5].0, 60.0);
205        assert!(segments[4].2.is_nan() && segments[5].2.is_nan());
206    }
207
208    #[test]
209    fn test_wind_sock_empty() {
210        let sock = WindSock::new(vec![]);
211        assert_eq!(sock.vector_for_range_stateless(50.0), Vector3::zeros());
212    }
213
214    #[test]
215    fn test_wind_sock_single_segment() {
216        // 16.0934 kmh (10 mph) @ 90° until 100m
217        let sock = WindSock::new(vec![(16.0934, 90.0, 100.0)]);
218
219        // Should have wind before 100m
220        let vec_50 = sock.vector_for_range_stateless(50.0);
221        println!("vec_50 = [{}, {}, {}]", vec_50[0], vec_50[1], vec_50[2]);
222        assert!(vec_50.norm() > 0.0);
223        // 90° wind from right (crosswind, McCoy): negative Z (lateral), zero Y, near-zero X (downrange)
224        assert!(
225            vec_50[2] < 0.0,
226            "Z (lateral) should be negative for 90° wind, got {}",
227            vec_50[2]
228        );
229        assert_eq!(vec_50[1], 0.0); // Zero Y component
230        assert!(
231            vec_50[0].abs() < 0.01,
232            "X (downrange) should be nearly zero for 90° wind, got {}",
233            vec_50[0]
234        );
235
236        // No wind after 100m
237        let vec_150 = sock.vector_for_range_stateless(150.0);
238        assert_eq!(vec_150, Vector3::zeros());
239    }
240
241    #[test]
242    fn test_wind_sock_multiple_segments() {
243        // Multiple wind segments (in kmh)
244        let sock = WindSock::new(vec![
245            (16.0934, 90.0, 50.0),  // 10 mph @ 90° until 50m
246            (24.1401, 45.0, 100.0), // 15 mph @ 45° until 100m
247            (8.0467, 180.0, 200.0), // 5 mph @ 180° until 200m
248        ]);
249
250        // Test each segment
251        let vec_25 = sock.vector_for_range_stateless(25.0);
252        println!("vec_25 = [{}, {}, {}]", vec_25[0], vec_25[1], vec_25[2]);
253        assert!(vec_25.norm() > 0.0);
254        assert!(vec_25[2] < 0.0, "90° wind should have negative Z (lateral)"); // 90° wind from right
255
256        let vec_75 = sock.vector_for_range_stateless(75.0);
257        println!("vec_75 = [{}, {}, {}]", vec_75[0], vec_75[1], vec_75[2]);
258        assert!(vec_75.norm() > vec_25.norm()); // 15 mph > 10 mph
259        assert!(vec_75[0] < 0.0); // 45° wind has negative X component
260        assert!(vec_75[2] < 0.0); // 45° wind has negative Z component
261
262        let vec_150 = sock.vector_for_range_stateless(150.0);
263        println!("vec_150 = [{}, {}, {}]", vec_150[0], vec_150[1], vec_150[2]);
264        assert!(vec_150.norm() < vec_75.norm()); // 5 mph < 15 mph
265        assert!(
266            vec_150[2].abs() < 0.01,
267            "180° wind should have near-zero Z (lateral), got {}",
268            vec_150[2]
269        ); // 180° wind (from behind)
270        assert!(
271            vec_150[0] > 0.0,
272            "180° wind should have positive X (tailwind, downrange), got {}",
273            vec_150[0]
274        ); // Tailwind
275
276        let vec_250 = sock.vector_for_range_stateless(250.0);
277        assert_eq!(vec_250, Vector3::zeros()); // Beyond all segments
278    }
279
280    #[test]
281    fn test_wind_conversion() {
282        // Test conversion: 16.0934 km/h = 4.47 m/s
283        let sock = WindSock::new(vec![(16.0934, 0.0, 100.0)]);
284        let vec = sock.vector_for_range_stateless(50.0);
285
286        let expected_speed = 16.0934 * KMH_TO_MPS;
287        assert!((vec.norm() - expected_speed).abs() < 0.01);
288    }
289
290    #[test]
291    fn test_wind_sock_boundary_is_upper_exclusive() {
292        // A segment's `until_distance_m` is exclusive: a query exactly at the
293        // boundary rolls to the next segment.
294        let sock = WindSock::new(vec![(16.0934, 90.0, 100.0), (32.1868, 270.0, 200.0)]);
295        // Just below 100 m -> first segment (90deg, negative Z).
296        assert!(sock.vector_for_range_stateless(99.999)[2] < 0.0);
297        // Exactly 100 m -> second segment (270deg, positive Z).
298        assert!(sock.vector_for_range_stateless(100.0)[2] > 0.0);
299        // Beyond the last boundary -> zero.
300        assert_eq!(sock.vector_for_range_stateless(200.0), Vector3::zeros());
301    }
302
303    #[test]
304    fn test_parse_wind_segment_str_units() {
305        // Imperial: 10 mph -> 16.0934 km/h, 100 yd -> 91.44 m.
306        let (kmh, ang, until) = parse_wind_segment_str("10:90:100", true).unwrap();
307        assert!((kmh - 16.09344).abs() < 1e-4);
308        assert_eq!(ang, 90.0);
309        assert!((until - 91.44).abs() < 1e-4);
310
311        // Metric: 5 m/s -> 18 km/h, 200 m stays 200 m.
312        let (kmh, ang, until) = parse_wind_segment_str("5:270:200", false).unwrap();
313        assert!((kmh - 18.0).abs() < 1e-9);
314        assert_eq!(ang, 270.0);
315        assert!((until - 200.0).abs() < 1e-9);
316
317        // Malformed inputs are rejected.
318        assert!(parse_wind_segment_str("10:90", true).is_err()); // too few fields
319        assert!(parse_wind_segment_str("10:bad:100", true).is_err()); // non-numeric
320        assert!(parse_wind_segment_str("10:90:0", true).is_err()); // zero until-distance
321        assert!(parse_wind_segment_str("-3:90:100", true).is_err()); // negative speed
322        // Non-finite values must be rejected (NaN comparisons would slip past < / <=).
323        assert!(parse_wind_segment_str("10:nan:5000", true).is_err());
324        assert!(parse_wind_segment_str("10:90:nan", true).is_err());
325        assert!(parse_wind_segment_str("inf:90:100", true).is_err());
326    }
327}