dive_deco/buhlmann/
compartment.rs

1use super::zhl_values::{ZHLParam, ZHLParams};
2use crate::{
3    common::{
4        powf, Depth, GradientFactor, InertGas, MbarPressure, PartialPressures, Pressure, RecordData,
5    },
6    BuhlmannConfig, Gas, Time,
7};
8#[cfg(feature = "serde")]
9use serde::{Deserialize, Serialize};
10
11#[derive(Copy, Clone, Debug, PartialEq)]
12#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
13pub struct Compartment {
14    // tissue number
15    pub no: u8,
16    // tolerable tissue ambient pressure
17    pub min_tolerable_amb_pressure: Pressure,
18    // helium saturation pressure
19    pub he_ip: Pressure,
20    // nitrogen saturation pressure
21    pub n2_ip: Pressure,
22    // total inert gas pressure (He + N2)
23    pub total_ip: Pressure,
24    // M-value (original)
25    pub m_value_raw: Pressure,
26    // M-value (calculated considering gradient factors)
27    pub m_value_calc: Pressure,
28    // compartment'a Buhlmann params (N2 half time, n2 'a' coefficient, n2 'b' coefficient, He half time, ..)
29    pub params: ZHLParams,
30    // Buhlmann model config (gradient factors, surface pressure)
31    model_config: BuhlmannConfig,
32}
33
34#[derive(Debug, PartialEq, Clone)]
35#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
36pub struct Supersaturation {
37    pub gf_99: f64,
38    pub gf_surf: f64,
39}
40
41impl Compartment {
42    pub fn new(no: u8, params: ZHLParams, model_config: BuhlmannConfig) -> Self {
43        let init_gas = Gas::air();
44        let init_gas_compound_pressures =
45            init_gas.inspired_partial_pressures(Depth::zero(), model_config.surface_pressure);
46        let n2_ip = init_gas_compound_pressures.n2;
47        let he_ip = init_gas_compound_pressures.he;
48
49        let mut compartment = Self {
50            no,
51            params,
52            n2_ip,
53            he_ip,
54            total_ip: he_ip + n2_ip,
55            m_value_raw: 0.,  // initial, recalculated later
56            m_value_calc: 0., // initial, recalculated later
57            min_tolerable_amb_pressure: 0.,
58            model_config,
59        };
60
61        // calculate initial minimal tolerable ambient pressure
62        let (_, gf_high) = compartment.model_config.gf;
63        compartment.m_value_raw = compartment.m_value(
64            Depth::zero(),
65            compartment.model_config.surface_pressure,
66            100,
67        );
68        compartment.m_value_calc = compartment.m_value_raw;
69        compartment.min_tolerable_amb_pressure = compartment.min_tolerable_amb_pressure(gf_high);
70
71        compartment
72    }
73
74    // recalculate tissue inert gasses saturation and tolerable pressure
75    pub fn recalculate(
76        &mut self,
77        record: &RecordData,
78        max_gf: GradientFactor,
79        surface_pressure: MbarPressure,
80    ) {
81        let (he_inert_pressure, n2_inert_pressure) =
82            self.compartment_inert_pressure(record, surface_pressure);
83
84        self.he_ip = he_inert_pressure;
85        self.n2_ip = n2_inert_pressure;
86        self.total_ip = he_inert_pressure + n2_inert_pressure;
87
88        // @todo m_value tuple
89        self.m_value_raw = self.m_value(record.depth, surface_pressure, 100);
90        self.m_value_calc = self.m_value(record.depth, surface_pressure, max_gf);
91
92        self.min_tolerable_amb_pressure = self.min_tolerable_amb_pressure(max_gf);
93    }
94
95    // tissue ceiling as depth
96    pub fn ceiling(&self) -> Depth {
97        let mut ceil = (self.min_tolerable_amb_pressure
98            - (self.model_config.surface_pressure as f64 / 1000.))
99            * 10.;
100        // cap ceiling at 0 if min tolerable leading compartment pressure depth equivalent negative
101        if ceil < 0. {
102            ceil = 0.;
103        }
104
105        Depth::from_meters(ceil)
106    }
107
108    // tissue supersaturation (gf99, surface gf)
109    pub fn supersaturation(&self, surface_pressure: MbarPressure, depth: Depth) -> Supersaturation {
110        let p_surf = (surface_pressure as f64) / 1000.;
111        let p_amb = p_surf + (depth.as_meters() / 10.);
112        let m_value = self.m_value_raw;
113        let m_value_surf = self.m_value(Depth::zero(), surface_pressure, 100);
114        let gf_99 = ((self.total_ip - p_amb) / (m_value - p_amb)) * 100.;
115        let gf_surf = ((self.total_ip - p_surf) / (m_value_surf - p_surf)) * 100.;
116
117        Supersaturation { gf_99, gf_surf }
118    }
119
120    fn m_value(
121        &self,
122        depth: Depth,
123        surface_pressure: MbarPressure,
124        max_gf: GradientFactor,
125    ) -> Pressure {
126        let weighted_zhl_params = self.weighted_zhl_params(self.he_ip, self.n2_ip);
127        let (_, a_coeff_adjusted, b_coeff_adjusted) =
128            self.max_gf_adjusted_zhl_params(weighted_zhl_params, max_gf);
129        let p_surf = (surface_pressure as f64) / 1000.;
130        let p_amb = p_surf + (depth.as_meters() / 10.);
131
132        a_coeff_adjusted + (p_amb / b_coeff_adjusted)
133    }
134
135    // tissue inert gasses pressure after record
136    fn compartment_inert_pressure(
137        &self,
138        record: &RecordData,
139        surface_pressure: MbarPressure,
140    ) -> (Pressure, Pressure) {
141        // (he, n2)
142        let RecordData { depth, time, gas } = record;
143        let PartialPressures {
144            n2: n2_pp,
145            he: he_pp,
146            ..
147        } = gas.inspired_partial_pressures(*depth, surface_pressure);
148
149        // partial pressure of inert gases in inspired gas (adjusted alveoli water vapor pressure)
150        let he_inspired_pp = he_pp;
151        let n2_inspired = n2_pp;
152
153        // tissue saturation pressure change for inert gasses
154        let (n2_half_time, _, _, he_half_time, ..) = self.params;
155        let he_p_comp_delta = self.compartment_pressure_delta_haldane(
156            InertGas::Helium,
157            he_inspired_pp,
158            *time,
159            he_half_time,
160        );
161        let n2_p_comp_delta = self.compartment_pressure_delta_haldane(
162            InertGas::Nitrogen,
163            n2_inspired,
164            *time,
165            n2_half_time,
166        );
167
168        // inert gasses pressures after applying delta P
169        let he_final = self.he_ip + he_p_comp_delta;
170        let n2_final = self.n2_ip + n2_p_comp_delta;
171
172        (he_final, n2_final)
173    }
174
175    // compartment pressure change for inert gas (Haldane equation)
176    fn compartment_pressure_delta_haldane(
177        &self,
178        inert_gas: InertGas,
179        gas_inspired_p: Pressure,
180        time: Time,
181        half_time: ZHLParam,
182    ) -> Pressure {
183        let inert_gas_load = match inert_gas {
184            InertGas::Helium => self.he_ip,
185            InertGas::Nitrogen => self.n2_ip,
186        };
187
188        // (Pi - Po)(1 - e^(-0.693t/half-time))
189        let factor = 1. - powf(2.0, -(time.as_minutes()) / half_time);
190
191        (gas_inspired_p - inert_gas_load) * factor
192    }
193
194    // tissue tolerable ambient pressure using GF slope, weighted Buhlmann ZHL params based on tissue inert gasses saturation proportions
195    fn min_tolerable_amb_pressure(&self, max_gf: GradientFactor) -> Pressure {
196        let weighted_zhl_params = self.weighted_zhl_params(self.he_ip, self.n2_ip);
197        let (_, a_coefficient_adjusted, b_coefficient_adjusted) =
198            self.max_gf_adjusted_zhl_params(weighted_zhl_params, max_gf);
199
200        (self.total_ip - a_coefficient_adjusted) * b_coefficient_adjusted
201    }
202
203    // weighted ZHL params (half time, a coefficient, b coefficient) based on N2 and He params and inert gasses proportions in tissue
204    pub fn weighted_zhl_params(
205        &self,
206        he_pp: Pressure,
207        n2_pp: Pressure,
208    ) -> (ZHLParam, ZHLParam, ZHLParam) {
209        fn weighted_param(
210            he_param: ZHLParam,
211            he_pp: Pressure,
212            n2_param: ZHLParam,
213            n2_pp: Pressure,
214        ) -> ZHLParam {
215            ((he_param * he_pp) + (n2_param * n2_pp)) / (he_pp + n2_pp)
216        }
217        let (n2_half_time, n2_a_coeff, n2_b_coeff, he_half_time, he_a_coeff, he_b_coeff) =
218            self.params;
219        (
220            weighted_param(he_half_time, he_pp, n2_half_time, n2_pp),
221            weighted_param(he_a_coeff, he_pp, n2_a_coeff, n2_pp),
222            weighted_param(he_b_coeff, he_pp, n2_b_coeff, n2_pp),
223        )
224    }
225
226    // adjust zhl params based on max gf
227    fn max_gf_adjusted_zhl_params(
228        &self,
229        params: (ZHLParam, ZHLParam, ZHLParam),
230        max_gf: GradientFactor,
231    ) -> (ZHLParam, ZHLParam, ZHLParam) {
232        let (half_time, a_coeff, b_coeff) = params;
233        let max_gf_fraction = max_gf as f64 / 100.;
234        let a_coefficient_adjusted = a_coeff * max_gf_fraction;
235        let b_coefficient_adjusted =
236            b_coeff / (max_gf_fraction - (max_gf_fraction * b_coeff) + b_coeff);
237
238        (half_time, a_coefficient_adjusted, b_coefficient_adjusted)
239    }
240}
241
242#[cfg(test)]
243mod tests {
244    use super::*;
245    use crate::{common::Gas, Time};
246
247    fn comp_1() -> Compartment {
248        let comp_1_params = (4., 1.2599, 0.5050, 1.51, 01.7424, 0.4245);
249        Compartment::new(1, comp_1_params, BuhlmannConfig::default())
250    }
251
252    fn comp_5() -> Compartment {
253        let comp_5_params = (27., 0.6200, 0.8126, 10.21, 0.9220, 0.7582);
254        Compartment::new(5, comp_5_params, BuhlmannConfig::default())
255    }
256
257    #[test]
258    fn test_constructor() {
259        let comp = comp_1();
260        assert_eq!(
261            comp,
262            Compartment {
263                no: 1,
264                min_tolerable_amb_pressure: -0.257127315,
265                he_ip: 0.0,
266                n2_ip: 0.750737,
267                total_ip: 0.750737,
268                m_value_raw: 3.265840594059406,
269                m_value_calc: 3.265840594059406,
270                params: (4.0, 1.2599, 0.505, 1.51, 1.7424, 0.4245),
271                // mocked config and state
272                model_config: BuhlmannConfig::default(),
273            }
274        );
275    }
276
277    #[test]
278    fn test_m_value_raw() {
279        let mut comp_1 = comp_1();
280        let mut comp_5 = comp_5();
281        let air = Gas::new(0.21, 0.);
282        let record = RecordData {
283            depth: Depth::zero(),
284            time: Time::from_seconds(1.),
285            gas: &air,
286        };
287        comp_1.recalculate(&record, 100, 1000);
288        comp_5.recalculate(&record, 100, 1000);
289        assert_eq!(comp_1.m_value_raw, 3.24009801980198);
290        assert_eq!(comp_5.m_value_raw, 1.8506177701206004);
291    }
292
293    #[test]
294    fn test_m_value_calc() {
295        let mut comp_1 = comp_1();
296        let mut comp_5 = comp_5();
297        let air = Gas::new(0.21, 0.);
298        let record = RecordData {
299            depth: Depth::zero(),
300            time: Time::from_seconds(1.),
301            gas: &air,
302        };
303        comp_1.recalculate(&record, 70, 1000);
304        comp_5.recalculate(&record, 70, 1000);
305        assert_eq!(comp_1.m_value_calc, 2.568068613861386);
306        assert_eq!(comp_5.m_value_calc, 1.5954324390844203);
307    }
308
309    #[test]
310    fn test_recalculation_ongassing() {
311        let mut comp = comp_5();
312        let air = Gas::new(0.21, 0.);
313        let record = RecordData {
314            depth: Depth::from_meters(30.),
315            time: Time::from_minutes(10.),
316            gas: &air,
317        };
318        comp.recalculate(&record, 100, 1000);
319        assert_eq!(comp.total_ip, 1.2850179204911072);
320    }
321
322    #[test]
323    fn test_weighted_params_trimix() {
324        let comp = comp_1();
325        let weighted_params = comp.weighted_zhl_params(0.5, 1. - (0.18 + 0.5));
326        assert_eq!(
327            weighted_params,
328            (2.481707317073171, 1.5541073170731705, 0.4559146341463414)
329        );
330    }
331
332    #[test]
333    fn test_min_pressure_calculation() {
334        let mut comp = comp_5();
335        let air = Gas::new(0.21, 0.);
336        let recprd = RecordData {
337            depth: Depth::from_meters(30.),
338            time: Time::from_minutes(10.),
339            gas: &air,
340        };
341        comp.recalculate(&recprd, 100, 100);
342        let min_tolerable_pressure = comp.min_tolerable_amb_pressure;
343        assert_eq!(min_tolerable_pressure, 0.40957969932131577);
344    }
345}
dive_deco/buhlmann/compartment.rs

dive_deco/buhlmann/
compartment.rs
