winds_engine/
field.rs

1//! The [`WindField`] trait and the concrete field types.
2
3use crate::{Error, Result};
4
5/// World-frame air velocity, m/s: +X is downrange (tailwind), +Y is up, +Z
6/// is the horizontal cross axis. Vertical wind is zero for the horizontal
7/// field types in this crate.
8#[derive(Debug, Clone, Copy, PartialEq)]
9pub struct Wind {
10    /// Downrange component, m/s.
11    pub x_mps: f64,
12    /// Vertical component, m/s.
13    pub y_mps: f64,
14    /// Cross-range component, m/s.
15    pub z_mps: f64,
16}
17
18impl Wind {
19    /// Calm air.
20    pub const ZERO: Self = Self {
21        x_mps: 0.0,
22        y_mps: 0.0,
23        z_mps: 0.0,
24    };
25
26    /// Horizontal wind from a speed and a meteorological direction: the
27    /// direction the air comes FROM, degrees clockwise, where 0° means from
28    /// downrange (+X, so the wind blows toward −X, a headwind) and 90° means
29    /// from +Z. Vertical component is zero.
30    pub fn from_speed_direction(speed_mps: f64, direction_from_deg: f64) -> Self {
31        let toward = (direction_from_deg + 180.0).to_radians();
32        Self {
33            x_mps: speed_mps * toward.cos(),
34            y_mps: 0.0,
35            z_mps: speed_mps * toward.sin(),
36        }
37    }
38
39    /// Horizontal speed.
40    pub fn speed_mps(&self) -> f64 {
41        self.x_mps.hypot(self.z_mps)
42    }
43
44    /// Meteorological direction the air comes FROM, degrees in [0, 360).
45    pub fn direction_from_deg(&self) -> f64 {
46        let toward = self.z_mps.atan2(self.x_mps).to_degrees() + 180.0;
47        toward.rem_euclid(360.0)
48    }
49}
50
51/// A deterministic, stateless wind field: same time and position always
52/// give the same wind, so integrator substeps, snapshots and seeded
53/// reproduction are exact.
54pub trait WindField {
55    /// World-frame air velocity at a time (seconds) and position
56    /// (`[x, y, z]` metres, Y up).
57    fn wind(&self, time_s: f64, position_m: [f64; 3]) -> Wind;
58}
59
60/// Uniform wind everywhere.
61#[derive(Debug, Clone, Copy)]
62pub struct ConstantField {
63    /// The constant air velocity.
64    pub wind: Wind,
65}
66
67impl WindField for ConstantField {
68    fn wind(&self, _time_s: f64, _position_m: [f64; 3]) -> Wind {
69        self.wind
70    }
71}
72
73/// One level of a vertical profile.
74#[derive(Debug, Clone, Copy, PartialEq)]
75pub struct ProfileLevel {
76    /// Height above the profile datum, metres.
77    pub altitude_m: f64,
78    /// Horizontal wind speed, m/s (nonnegative).
79    pub speed_mps: f64,
80    /// Meteorological direction the air comes FROM, degrees (see
81    /// [`Wind::from_speed_direction`]).
82    pub direction_from_deg: f64,
83}
84
85/// Piecewise-linear vertical profile: speed interpolates linearly in
86/// altitude, direction along the shortest circular arc. Below the first
87/// level and above the last, the nearest level holds.
88#[derive(Debug, Clone)]
89pub struct ProfileField {
90    levels: Vec<ProfileLevel>,
91}
92
93impl ProfileField {
94    /// Build from levels, which must be finite, nonnegative in speed and
95    /// strictly ascending in altitude.
96    pub fn new(mut levels: Vec<ProfileLevel>) -> Result<Self> {
97        if levels.is_empty() {
98            return Err(Error::invalid("profile", "needs at least one level"));
99        }
100        levels.sort_by(|a, b| a.altitude_m.total_cmp(&b.altitude_m));
101        for pair in levels.windows(2) {
102            if pair[0].altitude_m == pair[1].altitude_m {
103                return Err(Error::invalid(
104                    "profile",
105                    "altitudes must be strictly ascending",
106                ));
107            }
108        }
109        for level in &levels {
110            if !level.altitude_m.is_finite()
111                || !level.speed_mps.is_finite()
112                || !level.direction_from_deg.is_finite()
113            {
114                return Err(Error::invalid("profile", "levels must be finite"));
115            }
116            if level.speed_mps < 0.0 {
117                return Err(Error::invalid("profile", "speeds must be nonnegative"));
118            }
119        }
120        Ok(Self { levels })
121    }
122
123    /// The sorted levels.
124    pub fn levels(&self) -> &[ProfileLevel] {
125        &self.levels
126    }
127
128    /// Wind at an altitude, with circular direction interpolation.
129    pub fn wind_at(&self, altitude_m: f64) -> Wind {
130        let levels = &self.levels;
131        if altitude_m <= levels[0].altitude_m {
132            let l = levels[0];
133            return Wind::from_speed_direction(l.speed_mps, l.direction_from_deg);
134        }
135        for pair in levels.windows(2) {
136            let (lo, hi) = (pair[0], pair[1]);
137            if altitude_m <= hi.altitude_m {
138                let f = (altitude_m - lo.altitude_m) / (hi.altitude_m - lo.altitude_m);
139                let speed = lo.speed_mps + f * (hi.speed_mps - lo.speed_mps);
140                let mut delta = (hi.direction_from_deg - lo.direction_from_deg) % 360.0;
141                if delta > 180.0 {
142                    delta -= 360.0;
143                } else if delta < -180.0 {
144                    delta += 360.0;
145                }
146                let direction = lo.direction_from_deg + f * delta;
147                return Wind::from_speed_direction(speed, direction);
148            }
149        }
150        let l = levels[levels.len() - 1];
151        Wind::from_speed_direction(l.speed_mps, l.direction_from_deg)
152    }
153}
154
155impl WindField for ProfileField {
156    fn wind(&self, _time_s: f64, position_m: [f64; 3]) -> Wind {
157        self.wind_at(position_m[1])
158    }
159}
winds_engine/field.rs

winds_engine/
field.rs
