hpr_sim/
environment.rs

1//! The flight's surroundings: the Earth model with the launch site, the atmosphere and the wind.
2
3use std::sync::Arc;
4
5use hpr_atmos::{AirState, Atmosphere, AtmosphereModel, ConstantWind, Wind};
6use hpr_core::DVec3;
7use hpr_core::earth::Earth;
8use hpr_core::geodesy::Geodetic;
9
10use crate::error::SimError;
11
12/// The Earth, atmosphere and wind a flight sees.
13///
14/// Heights: the state and the Earth model use the launch frame and ellipsoidal heights. The
15/// atmosphere and wind take height above mean sea level, `H = h − N`, with the geoid undulation
16/// `N` at the site given here (hpr has no geoid model; `docs/physics/geodesy.md`). The ground is
17/// the ellipsoidal height of the site. Wind vectors are taken in the launch frame's axes.
18#[derive(Debug, Clone)]
19pub struct Environment {
20    /// The Earth: the launch site and frame, gravity and rotation.
21    pub earth: Earth,
22    /// The geoid undulation `N` at the site, m (`h = H + N`).
23    pub geoid_undulation_m: f64,
24    /// The atmosphere, by height above mean sea level. Shared, so environments clone cheaply
25    /// across threads.
26    pub atmosphere: Arc<dyn Atmosphere>,
27    /// The wind, by height above mean sea level, in launch-frame axes.
28    pub wind: Arc<dyn Wind>,
29}
30
31impl Environment {
32    /// An environment from its parts, with `N = 0`.
33    pub fn new(
34        earth: Earth,
35        atmosphere: impl Atmosphere + 'static,
36        wind: impl Wind + 'static,
37    ) -> Self {
38        Self {
39            earth,
40            geoid_undulation_m: 0.0,
41            atmosphere: Arc::new(atmosphere),
42            wind: Arc::new(wind),
43        }
44    }
45
46    /// WGS 84 at `site` (ellipsoidal height), the 1976 standard atmosphere and no wind.
47    ///
48    /// # Errors
49    ///
50    /// [`SimError::Core`] for an invalid site.
51    pub fn standard(site: Geodetic) -> Result<Self, SimError> {
52        Ok(Self::new(
53            Earth::wgs84(site)?,
54            AtmosphereModel::default(),
55            ConstantWind::calm(),
56        ))
57    }
58
59    /// The same environment with geoid undulation `N`, m.
60    #[must_use]
61    pub fn with_geoid_undulation_m(mut self, undulation_m: f64) -> Self {
62        self.geoid_undulation_m = undulation_m;
63        self
64    }
65
66    /// The same environment with `wind` in place of its wind.
67    #[must_use]
68    pub fn with_wind(mut self, wind: impl Wind + 'static) -> Self {
69        self.wind = Arc::new(wind);
70        self
71    }
72
73    /// The launch site.
74    #[must_use]
75    pub fn site(&self) -> Geodetic {
76        self.earth.frame().origin()
77    }
78
79    /// The wind's velocity at `height_msl_m` above mean sea level, in the launch frame's
80    /// East-North-Up axes, m/s. Every place the flight reads the wind reads it here, so a wind of
81    /// a program's own that returns a velocity that is not finite (NaN or infinite) is refused
82    /// where it is read, naming the wind and the height, rather than surfacing later as some other
83    /// quantity that isn't finite (issue #237).
84    ///
85    /// # Errors
86    ///
87    /// [`SimError::Atmosphere`] if the wind refuses the height; [`SimError::Domain`] with
88    /// [`WIND_NOT_FINITE`] and the height above mean sea level, m, if any component of the
89    /// velocity is not finite.
90    pub(crate) fn wind_enu_m_s(&self, height_msl_m: f64) -> Result<DVec3, SimError> {
91        let velocity_enu_m_s = self.wind.wind(height_msl_m)?.velocity_enu_m_s;
92        if velocity_enu_m_s.is_finite() {
93            Ok(velocity_enu_m_s)
94        } else {
95            Err(SimError::Domain {
96                what: WIND_NOT_FINITE,
97                value: height_msl_m,
98            })
99        }
100    }
101
102    /// The air at `height_msl_m` above mean sea level. Every place the flight reads the air reads
103    /// it here, so an atmosphere of a program's own that returns air the flight can't use is
104    /// refused where it is read, naming the field and the height, rather than flying on with a
105    /// wrong number (issue #301). Before, a density that was NaN or negative turned the drag off
106    /// without an error: a separated body landed at about 140 m/s, a climb went twice as high.
107    ///
108    /// Every field must be finite. The temperature, speed of sound and viscosity must also be
109    /// positive. The density and pressure may be zero, as in a vacuum or far above the 1976
110    /// standard atmosphere's top, where both shrink toward zero, but not negative. The air is
111    /// returned unchanged, so good air flies exactly as it did before the check.
112    ///
113    /// # Errors
114    ///
115    /// [`SimError::Atmosphere`] if the atmosphere refuses the height; [`SimError::Domain`] with
116    /// the field's constant ([`AIR_DENSITY_REFUSED`], [`AIR_PRESSURE_REFUSED`],
117    /// [`AIR_TEMPERATURE_REFUSED`], [`AIR_SPEED_OF_SOUND_REFUSED`] or
118    /// [`AIR_VISCOSITY_REFUSED`]) and the height above mean sea level, m, for the first field,
119    /// in that order, the flight can't use.
120    pub(crate) fn air_at(&self, height_msl_m: f64) -> Result<AirState, SimError> {
121        let air = self.atmosphere.air(height_msl_m)?.air;
122        let not_negative = |x: f64| x.is_finite() && x >= 0.0;
123        let positive = |x: f64| x.is_finite() && x > 0.0;
124        let refused = if !not_negative(air.density_kg_m3) {
125            AIR_DENSITY_REFUSED
126        } else if !not_negative(air.pressure_pa) {
127            AIR_PRESSURE_REFUSED
128        } else if !positive(air.temperature_k) {
129            AIR_TEMPERATURE_REFUSED
130        } else if !positive(air.speed_of_sound_m_s) {
131            AIR_SPEED_OF_SOUND_REFUSED
132        } else if !positive(air.dynamic_viscosity_pa_s) {
133            AIR_VISCOSITY_REFUSED
134        } else {
135            return Ok(air);
136        };
137        Err(SimError::Domain {
138            what: refused,
139            value: height_msl_m,
140        })
141    }
142}
143
144/// What a [`SimError::Domain`] names when the wind's velocity is not finite; its value is the
145/// height above mean sea level, m, at which the wind was read.
146pub(crate) const WIND_NOT_FINITE: &str =
147    "height above sea level, m, at which the wind's velocity is not finite";
148
149/// What a [`SimError::Domain`] names when the air's density is negative or not finite; its value
150/// is the height above mean sea level, m, at which the air was read.
151pub(crate) const AIR_DENSITY_REFUSED: &str =
152    "height above sea level, m, at which the air's density is negative or not finite";
153
154/// What a [`SimError::Domain`] names when the air's pressure is negative or not finite; its value
155/// is the height above mean sea level, m, at which the air was read.
156pub(crate) const AIR_PRESSURE_REFUSED: &str =
157    "height above sea level, m, at which the air's pressure is negative or not finite";
158
159/// What a [`SimError::Domain`] names when the air's temperature is zero, negative or not finite;
160/// its value is the height above mean sea level, m, at which the air was read.
161pub(crate) const AIR_TEMPERATURE_REFUSED: &str =
162    "height above sea level, m, at which the air's temperature is zero, negative or not finite";
163
164/// What a [`SimError::Domain`] names when the air's speed of sound is zero, negative or not
165/// finite; its value is the height above mean sea level, m, at which the air was read.
166pub(crate) const AIR_SPEED_OF_SOUND_REFUSED: &str = "height above sea level, m, at which the \
167     air's speed of sound is zero, negative or not finite";
168
169/// What a [`SimError::Domain`] names when the air's dynamic viscosity is zero, negative or not
170/// finite; its value is the height above mean sea level, m, at which the air was read.
171pub(crate) const AIR_VISCOSITY_REFUSED: &str = "height above sea level, m, at which the air's \
172     viscosity is zero, negative or not finite";
173
174#[cfg(test)]
175mod tests {
176    use hpr_atmos::ConstantWind;
177
178    use std::sync::Arc;
179
180    use hpr_atmos::Atmosphere;
181
182    use super::Environment;
183    use crate::testing::site;
184
185    #[test]
186    fn with_wind_replaces_only_the_wind() {
187        let calm = Environment::standard(site()).unwrap();
188        // 5 m/s from the west blows toward the east.
189        let windy = calm
190            .clone()
191            .with_wind(ConstantWind::new(5.0, 1.5 * std::f64::consts::PI).unwrap());
192        let at =
193            |environment: &Environment| environment.wind.wind(1500.0).unwrap().velocity_enu_m_s;
194        assert_eq!(at(&calm).length(), 0.0);
195        assert!((at(&windy).x - 5.0).abs() < 1e-12 && at(&windy).y.abs() < 1e-12);
196        assert_eq!(windy.site(), calm.site());
197        assert_eq!(windy.geoid_undulation_m, calm.geoid_undulation_m);
198        assert_eq!(
199            windy.atmosphere.air(1500.0).unwrap(),
200            calm.atmosphere.air(1500.0).unwrap()
201        );
202    }
203
204    /// The check on the air refuses nothing hpr's own atmospheres return, from 5 km below sea
205    /// level to a million kilometers up, where the standard's pressure and density have shrunk to
206    /// zero: the standard, offset as far as it allows either way, a humid sounding extrapolated
207    /// past both ends, and the tests' vacuum. It returns the air unchanged.
208    #[test]
209    fn the_check_on_the_air_takes_every_stock_atmosphere() {
210        use hpr_atmos::{AtmosphereModel, Ussa76};
211
212        use crate::testing::UniformAir;
213
214        let sounding: AtmosphereModel = serde_json::from_str(
215            r#"{"model":"sounding","latitude_rad":0.6,"levels":[
216            {"height_msl_m":1400.0,"temperature_k":300.0,"pressure_pa":85000.0,
217             "relative_humidity":1.0},
218            {"height_msl_m":3000.0,"temperature_k":290.0,"pressure_pa":70000.0,
219             "relative_humidity":0.5}]}"#,
220        )
221        .unwrap();
222        let atmospheres: [Arc<dyn Atmosphere>; 5] = [
223            Arc::new(Ussa76::standard()),
224            Arc::new(Ussa76::with_offset(-186.0, 30_000.0).unwrap()),
225            Arc::new(Ussa76::with_offset(60.0, 110_000.0).unwrap()),
226            Arc::new(sounding),
227            Arc::new(UniformAir::vacuum()),
228        ];
229        // Every 500 m from 5 km below sea level to 1 km above it, then 25% higher each time, and
230        // last a million kilometers up.
231        let mut heights_msl_m: Vec<f64> = (-10..=2).map(|i| f64::from(i) * 500.0).collect();
232        while let Some(&last) = heights_msl_m.last().filter(|&&h| h < 1.0e9) {
233            heights_msl_m.push((last * 1.25).min(1.0e9));
234        }
235        assert_eq!(heights_msl_m.last(), Some(&1.0e9));
236        for atmosphere in atmospheres {
237            let mut environment = Environment::standard(site()).unwrap();
238            environment.atmosphere = atmosphere;
239            for &height_msl_m in &heights_msl_m {
240                let air = environment.air_at(height_msl_m).unwrap();
241                assert_eq!(
242                    air,
243                    environment.atmosphere.air(height_msl_m).unwrap().air,
244                    "{height_msl_m}"
245                );
246            }
247        }
248        // The standard's own air at the top of the sweep has a density and pressure of zero, so
249        // the check must take a zero.
250        let far = Environment::standard(site())
251            .unwrap()
252            .air_at(1.0e9)
253            .unwrap();
254        assert_eq!((far.density_kg_m3, far.pressure_pa), (0.0, 0.0));
255    }
256
257    /// Each field the flight can't use is refused by name, with the height; the first such field,
258    /// in the order density, pressure, temperature, speed of sound, viscosity, is the one named.
259    #[test]
260    fn the_check_on_the_air_names_the_field_and_the_height() {
261        use super::{
262            AIR_DENSITY_REFUSED, AIR_PRESSURE_REFUSED, AIR_SPEED_OF_SOUND_REFUSED,
263            AIR_TEMPERATURE_REFUSED, AIR_VISCOSITY_REFUSED,
264        };
265        use crate::error::SimError;
266        use crate::testing::UniformAir;
267
268        let refusals = [
269            (AIR_DENSITY_REFUSED, true),
270            (AIR_PRESSURE_REFUSED, true),
271            (AIR_TEMPERATURE_REFUSED, false),
272            (AIR_SPEED_OF_SOUND_REFUSED, false),
273            (AIR_VISCOSITY_REFUSED, false),
274        ];
275        let set = |air: &mut hpr_atmos::AirState, field: usize, value: f64| {
276            *[
277                &mut air.density_kg_m3,
278                &mut air.pressure_pa,
279                &mut air.temperature_k,
280                &mut air.speed_of_sound_m_s,
281                &mut air.dynamic_viscosity_pa_s,
282            ]
283            .into_iter()
284            .nth(field)
285            .unwrap() = value;
286        };
287        let check = |air: hpr_atmos::AirState| {
288            let mut environment = Environment::standard(site()).unwrap();
289            environment.atmosphere = Arc::new(UniformAir(air));
290            environment.air_at(1234.5)
291        };
292        for (field, (refusal, takes_zero)) in refusals.into_iter().enumerate() {
293            let mut bad = vec![
294                f64::NAN,
295                f64::INFINITY,
296                f64::NEG_INFINITY,
297                -1.0,
298                -f64::MIN_POSITIVE,
299            ];
300            if !takes_zero {
301                bad.extend([0.0, -0.0]);
302            }
303            for value in bad {
304                // Every later field spoilt too: the first is the one named.
305                let mut air = UniformAir::sea_level().0;
306                for later in field..5 {
307                    set(&mut air, later, value);
308                }
309                match check(air) {
310                    Err(SimError::Domain {
311                        what,
312                        value: height,
313                    }) => {
314                        assert_eq!(what, refusal, "{field} {value}");
315                        assert_eq!(height, 1234.5, "{field} {value}");
316                    }
317                    other => panic!("{field} {value}: {other:?}"),
318                }
319            }
320            let mut tiny = UniformAir::sea_level().0;
321            set(&mut tiny, field, f64::MIN_POSITIVE);
322            assert_eq!(check(tiny).unwrap(), tiny, "{field}");
323            if takes_zero {
324                for zero in [0.0, -0.0] {
325                    let mut air = UniformAir::sea_level().0;
326                    set(&mut air, field, zero);
327                    assert_eq!(check(air).unwrap(), air, "{field} {zero}");
328                }
329            }
330        }
331    }
332}
hpr_sim/environment.rs

hpr_sim/
environment.rs
