zhl16/
body_model.rs

1use crate::gas::{Gas, GasMixture};
2use crate::quantities::{Frequency, FrequencyUnit, Pressure, PressureUnit, Time, TimeUnit};
3use crate::tissue_constants::TISSUE_CONSTANTS;
4use core::array::from_fn;
5use core::fmt;
6
7const WATER_VAPOR_PRESSURE: Pressure = Pressure::new_unchecked(0.0627, PressureUnit::Bar);
8
9#[derive(Copy, Clone, Debug, Eq, PartialEq)]
10pub enum BodyModelError {
11    NonPositiveTimeStep,
12    InvalidGradientFactor,
13}
14
15impl fmt::Display for BodyModelError {
16    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
17        match self {
18            BodyModelError::NonPositiveTimeStep => f.write_str("time step must be positive"),
19            BodyModelError::InvalidGradientFactor => {
20                f.write_str("gradient factor must be between 0.0 and 1.0")
21            }
22        }
23    }
24}
25
26#[derive(Copy, Clone)]
27pub(crate) struct TissueModel {
28    gases: [TissueModelPerGas; 2],
29}
30
31impl TissueModel {
32    pub(crate) const fn new(nitrogen: TissueModelPerGas, helium: TissueModelPerGas) -> Self {
33        Self {
34            gases: [nitrogen, helium],
35        }
36    }
37
38    #[cfg(test)]
39    pub(crate) fn get_gas(&self, gas: Gas) -> TissueModelPerGas {
40        self.gases[gas as usize]
41    }
42}
43
44#[derive(Copy, Clone)]
45pub(crate) struct TissueModelPerGas {
46    k: Frequency,
47    a: Pressure,
48    b: f64,
49}
50
51impl TissueModelPerGas {
52    pub const fn new(half_time: Time, a: Pressure, b: f64) -> Self {
53        Self {
54            k: Frequency::new_unchecked(
55                core::f64::consts::LN_2 / half_time.to(TimeUnit::Seconds),
56                FrequencyUnit::Hertz,
57            ),
58            a,
59            b,
60        }
61    }
62    #[cfg(test)]
63    pub(crate) fn get_k(&self) -> Frequency {
64        self.k
65    }
66    pub fn get_a(&self) -> Pressure {
67        self.a
68    }
69    pub fn get_b(&self) -> f64 {
70        self.b
71    }
72}
73
74#[derive(Copy, Clone)]
75pub struct CompartmentState {
76    pressures: [Pressure; 2],
77    tissue_model: TissueModel,
78}
79
80impl CompartmentState {
81    pub(crate) fn new(
82        nitrogen_pressures: Pressure,
83        helium_pressures: Pressure,
84        tissue_model: TissueModel,
85    ) -> Self {
86        let mut pressures = [Pressure::zero(); 2];
87        pressures[Gas::Nitrogen as usize] = nitrogen_pressures;
88        pressures[Gas::Helium as usize] = helium_pressures;
89        Self {
90            pressures,
91            tissue_model,
92        }
93    }
94
95    pub(crate) fn evolve(
96        &mut self,
97        constants: &TissueModel,
98        ambient_pressure: Pressure,
99        gas_fractions: &GasMixture,
100        dt: Time,
101    ) {
102        for gas in Gas::ALL {
103            let idx = gas as usize;
104            let tissue = &constants.gases[idx];
105            let pressure = &mut self.pressures[idx];
106            let fraction = gas_fractions.inert_fraction(gas);
107
108            let inspired = (ambient_pressure - WATER_VAPOR_PRESSURE) * fraction;
109            let rate = (inspired - *pressure) * tissue.k;
110            *pressure += rate * dt;
111        }
112    }
113
114    pub fn get_nitrogen_pressure(&self) -> Pressure {
115        self.pressures[Gas::Nitrogen as usize]
116    }
117
118    pub fn get_helium_pressure(&self) -> Pressure {
119        self.pressures[Gas::Helium as usize]
120    }
121
122    pub fn safe_ambient_pressure(&self, gf: f64) -> Result<Pressure, BodyModelError> {
123        if !(0.0..=1.0).contains(&gf) {
124            return Err(BodyModelError::InvalidGradientFactor);
125        }
126
127        let p_n2 = self.pressures[0];
128        let p_he = self.pressures[1];
129        let p_total = p_n2 + p_he;
130
131        if p_total < Pressure::new_unchecked(0.001, PressureUnit::Bar) {
132            return Ok(Pressure::zero());
133        }
134
135        let model_n2 = self.tissue_model.gases[0];
136        let model_he = self.tissue_model.gases[1];
137
138        let a = (p_n2 / p_total) * model_n2.get_a() + (p_he / p_total) * model_he.get_a();
139
140        let b = (model_n2.get_b() * p_n2 + model_he.get_b() * p_he) / p_total;
141
142        let numerator = p_total - a * gf;
143        let denominator = gf / b + (1.0 - gf);
144
145        Ok(numerator / denominator)
146    }
147}
148
149pub struct BodyState {
150    pub(crate) compartments: [CompartmentState; 16],
151}
152
153impl BodyState {
154    pub fn new(nitrogen_pressure: Pressure, helium_pressure: Pressure) -> Self {
155        Self {
156            compartments: from_fn(|i| {
157                CompartmentState::new(nitrogen_pressure, helium_pressure, TISSUE_CONSTANTS[i])
158            }),
159        }
160    }
161
162    pub fn compartments(&self) -> &[CompartmentState; 16] {
163        &self.compartments
164    }
165
166    pub fn evolve(
167        &mut self,
168        ambient_pressure: Pressure,
169        gas_mixture: &GasMixture,
170        dt: Time,
171    ) -> Result<(), BodyModelError> {
172        if !(dt.to(TimeUnit::Seconds) > 0.0) {
173            return Err(BodyModelError::NonPositiveTimeStep);
174        }
175
176        for (i, compartment) in self.compartments.iter_mut().enumerate() {
177            let constants = &TISSUE_CONSTANTS[i];
178            compartment.evolve(constants, ambient_pressure, gas_mixture, dt);
179        }
180
181        Ok(())
182    }
183}
zhl16/body_model.rs

zhl16/
body_model.rs
