zhl16 0.1.0

A no_std Rust implementation of core Bühlmann ZHL-16 tissue loading primitives.
Documentation
#[cfg(test)]
mod tests {
    use crate::body_model::*;
    use crate::gas::{Gas, GasMixture};
    use crate::quantities::*;
    use crate::tissue_constants::TISSUE_CONSTANTS;

    fn assert_close(actual: f64, expected: f64, epsilon: f64) {
        let diff = if actual > expected {
            actual - expected
        } else {
            expected - actual
        };

        assert!(
            diff <= epsilon,
            "expected {actual} to be within {epsilon} of {expected}"
        );
    }

    fn assert_pressure_close(actual: Pressure, expected_bar: f64, epsilon: f64) {
        assert_close(actual.to(PressureUnit::Bar), expected_bar, epsilon);
    }

    fn gas_mixture(oxygen_fraction: f64, helium_fraction: f64) -> GasMixture {
        GasMixture::new(oxygen_fraction, helium_fraction).expect("valid gas mixture")
    }

    fn expected_evolved_pressure(
        initial_bar: f64,
        ambient_bar: f64,
        fraction: f64,
        tissue: TissueModelPerGas,
        dt_seconds: f64,
    ) -> f64 {
        let inspired = (ambient_bar - 0.0627) * fraction;
        initial_bar
            + (inspired - initial_bar) * tissue.get_k().to(FrequencyUnit::Hertz) * dt_seconds
    }

    fn expected_safe_ambient_pressure(
        model: &TissueModel,
        nitrogen_pressure_bar: f64,
        helium_pressure_bar: f64,
        gf: f64,
    ) -> f64 {
        let total_pressure_bar = nitrogen_pressure_bar + helium_pressure_bar;

        if total_pressure_bar < 0.001 {
            return 0.0;
        }

        let nitrogen_model = model.get_gas(Gas::Nitrogen);
        let helium_model = model.get_gas(Gas::Helium);
        let a = (nitrogen_pressure_bar / total_pressure_bar)
            * nitrogen_model.get_a().to(PressureUnit::Bar)
            + (helium_pressure_bar / total_pressure_bar)
                * helium_model.get_a().to(PressureUnit::Bar);
        let b = (nitrogen_model.get_b() * nitrogen_pressure_bar
            + helium_model.get_b() * helium_pressure_bar)
            / total_pressure_bar;

        (total_pressure_bar - a * gf) / (gf / b + (1.0 - gf))
    }

    #[test]
    fn tissue_model_per_gas_stores_half_time_as_frequency() {
        let tissue = TissueModelPerGas::new(
            Time::new_unchecked(5.0, TimeUnit::Minutes),
            Pressure::new_unchecked(1.1696, PressureUnit::Bar),
            0.5578,
        );

        assert_close(
            tissue.get_k().to(FrequencyUnit::Hertz),
            core::f64::consts::LN_2 / 300.0,
            1E-12,
        );
        assert_eq!(tissue.get_a().to(PressureUnit::Bar), 1.1696);
        assert_eq!(tissue.get_b(), 0.5578);
    }

    #[test]
    fn body_state_initializes_all_compartments_with_same_pressures() {
        let body = BodyState::new(
            Pressure::new_unchecked(0.79, PressureUnit::Bar),
            Pressure::new_unchecked(0.0, PressureUnit::Bar),
        );

        for compartment in body.compartments {
            assert_eq!(
                compartment.get_nitrogen_pressure().to(PressureUnit::Bar),
                0.79
            );
            assert_eq!(compartment.get_helium_pressure().to(PressureUnit::Bar), 0.0);
        }
    }

    #[test]
    fn compartment_evolve_updates_each_gas_pressure() {
        let constants = TissueModel::new(
            TissueModelPerGas::new(
                Time::new_unchecked(10.0, TimeUnit::Minutes),
                Pressure::new_unchecked(1.0, PressureUnit::Bar),
                0.5,
            ),
            TissueModelPerGas::new(
                Time::new_unchecked(5.0, TimeUnit::Minutes),
                Pressure::new_unchecked(1.0, PressureUnit::Bar),
                0.5,
            ),
        );
        let mut compartment = CompartmentState::new(
            Pressure::new_unchecked(0.70, PressureUnit::Bar),
            Pressure::new_unchecked(0.10, PressureUnit::Bar),
            TISSUE_CONSTANTS[0],
        );

        compartment.evolve(
            &constants,
            Pressure::new_unchecked(2.0, PressureUnit::Bar),
            &gas_mixture(0.10, 0.20),
            Time::new_unchecked(60.0, TimeUnit::Seconds),
        );

        let inspired_n2 = (2.0 - 0.0627) * 0.70;
        let expected_n2 = 0.70 + (inspired_n2 - 0.70) * (core::f64::consts::LN_2 / 600.0) * 60.0;
        let inspired_he = (2.0 - 0.0627) * 0.20;
        let expected_he = 0.10 + (inspired_he - 0.10) * (core::f64::consts::LN_2 / 300.0) * 60.0;

        assert_close(
            compartment.get_nitrogen_pressure().to(PressureUnit::Bar),
            expected_n2,
            1E-12,
        );
        assert_close(
            compartment.get_helium_pressure().to(PressureUnit::Bar),
            expected_he,
            1E-12,
        );
    }

    #[test]
    fn compartment_evolve_with_zero_time_keeps_pressures() {
        let mut compartment = CompartmentState::new(
            Pressure::new_unchecked(0.70, PressureUnit::Bar),
            Pressure::new_unchecked(0.10, PressureUnit::Bar),
            TISSUE_CONSTANTS[0],
        );

        compartment.evolve(
            &TISSUE_CONSTANTS[0],
            Pressure::new_unchecked(3.0, PressureUnit::Bar),
            &gas_mixture(0.20, 0.30),
            Time::zero(),
        );

        assert_pressure_close(compartment.get_nitrogen_pressure(), 0.70, 1E-12);
        assert_pressure_close(compartment.get_helium_pressure(), 0.10, 1E-12);
    }

    #[test]
    fn body_evolve_uses_tissue_constants_for_each_compartment() {
        let mut body = BodyState::new(
            Pressure::new_unchecked(0.79, PressureUnit::Bar),
            Pressure::new_unchecked(0.0, PressureUnit::Bar),
        );

        body.evolve(
            Pressure::new_unchecked(2.0, PressureUnit::Bar),
            &gas_mixture(0.10, 0.20),
            Time::new_unchecked(60.0, TimeUnit::Seconds),
        )
        .expect("positive time step");

        let first = body.compartments[0]
            .get_nitrogen_pressure()
            .to(PressureUnit::Bar);
        let last = body.compartments[15]
            .get_nitrogen_pressure()
            .to(PressureUnit::Bar);
        let inspired_n2 = (2.0 - 0.0627) * 0.70;
        let expected_first = 0.79 + (inspired_n2 - 0.79) * (core::f64::consts::LN_2 / 300.0) * 60.0;

        assert_close(first, expected_first, 1E-12);
        assert!(
            first > last,
            "faster first compartment should move farther than the last compartment"
        );
    }

    #[test]
    fn body_evolve_updates_every_compartment_for_both_gases() {
        let mut body = BodyState::new(
            Pressure::new_unchecked(0.79, PressureUnit::Bar),
            Pressure::new_unchecked(0.12, PressureUnit::Bar),
        );

        body.evolve(
            Pressure::new_unchecked(3.4, PressureUnit::Bar),
            &gas_mixture(0.20, 0.30),
            Time::new_unchecked(90.0, TimeUnit::Seconds),
        )
        .expect("positive time step");

        for (i, model) in TISSUE_CONSTANTS.iter().enumerate() {
            let compartment = &body.compartments[i];
            let expected_nitrogen =
                expected_evolved_pressure(0.79, 3.4, 0.50, model.get_gas(Gas::Nitrogen), 90.0);
            let expected_helium =
                expected_evolved_pressure(0.12, 3.4, 0.30, model.get_gas(Gas::Helium), 90.0);

            assert_pressure_close(
                compartment.get_nitrogen_pressure(),
                expected_nitrogen,
                1E-12,
            );
            assert_pressure_close(compartment.get_helium_pressure(), expected_helium, 1E-12);
        }
    }

    #[test]
    fn body_evolve_rejects_non_positive_time_steps() {
        let mut body = BodyState::new(
            Pressure::new_unchecked(0.79, PressureUnit::Bar),
            Pressure::zero(),
        );

        assert_eq!(
            body.evolve(
                Pressure::new_unchecked(2.0, PressureUnit::Bar),
                &gas_mixture(0.10, 0.20),
                Time::zero(),
            ),
            Err(BodyModelError::NonPositiveTimeStep)
        );
        assert_eq!(
            body.evolve(
                Pressure::new_unchecked(2.0, PressureUnit::Bar),
                &gas_mixture(0.10, 0.20),
                Time::new_unchecked(-1.0, TimeUnit::Seconds),
            ),
            Err(BodyModelError::NonPositiveTimeStep)
        );
        assert_eq!(
            body.evolve(
                Pressure::new_unchecked(2.0, PressureUnit::Bar),
                &gas_mixture(0.10, 0.20),
                Time::new_unchecked(f64::NAN, TimeUnit::Seconds),
            ),
            Err(BodyModelError::NonPositiveTimeStep)
        );
    }

    #[test]
    fn body_state_assigns_matching_tissue_model_to_each_compartment() {
        let body = BodyState::new(
            Pressure::new_unchecked(1.20, PressureUnit::Bar),
            Pressure::new_unchecked(0.50, PressureUnit::Bar),
        );

        for (i, model) in TISSUE_CONSTANTS.iter().enumerate() {
            let expected = expected_safe_ambient_pressure(model, 1.20, 0.50, 0.85);

            assert_pressure_close(
                body.compartments[i]
                    .safe_ambient_pressure(0.85)
                    .expect("valid gradient factor"),
                expected,
                1E-12,
            );
        }
    }

    #[test]
    fn safe_ambient_pressure_handles_gradient_factor_bounds() {
        let compartment = CompartmentState::new(
            Pressure::new_unchecked(1.80, PressureUnit::Bar),
            Pressure::zero(),
            TISSUE_CONSTANTS[0],
        );
        let nitrogen_model = TISSUE_CONSTANTS[0].get_gas(Gas::Nitrogen);
        let expected_gf_one =
            (1.80 - nitrogen_model.get_a().to(PressureUnit::Bar)) * nitrogen_model.get_b();

        assert_pressure_close(
            compartment
                .safe_ambient_pressure(0.0)
                .expect("valid gradient factor"),
            1.80,
            1E-12,
        );
        assert_pressure_close(
            compartment
                .safe_ambient_pressure(1.0)
                .expect("valid gradient factor"),
            expected_gf_one,
            1E-12,
        );
    }

    #[test]
    fn safe_ambient_pressure_handles_mixed_gas_weighting() {
        let compartment = CompartmentState::new(
            Pressure::new_unchecked(1.20, PressureUnit::Bar),
            Pressure::new_unchecked(0.50, PressureUnit::Bar),
            TISSUE_CONSTANTS[0],
        );
        let expected = expected_safe_ambient_pressure(&TISSUE_CONSTANTS[0], 1.20, 0.50, 0.85);

        assert_pressure_close(
            compartment
                .safe_ambient_pressure(0.85)
                .expect("valid gradient factor"),
            expected,
            1E-12,
        );
    }

    #[test]
    fn safe_ambient_pressure_returns_zero_for_near_empty_compartment() {
        let compartment = CompartmentState::new(
            Pressure::new_unchecked(0.0004, PressureUnit::Bar),
            Pressure::new_unchecked(0.0005, PressureUnit::Bar),
            TISSUE_CONSTANTS[0],
        );

        assert_eq!(
            compartment
                .safe_ambient_pressure(0.85)
                .expect("valid gradient factor")
                .to(PressureUnit::Bar),
            0.0
        );
    }

    #[test]
    fn safe_ambient_pressure_rejects_invalid_gradient_factors() {
        let compartment = CompartmentState::new(
            Pressure::new_unchecked(1.20, PressureUnit::Bar),
            Pressure::new_unchecked(0.50, PressureUnit::Bar),
            TISSUE_CONSTANTS[0],
        );

        assert_eq!(
            compartment.safe_ambient_pressure(-0.01),
            Err(BodyModelError::InvalidGradientFactor)
        );
        assert_eq!(
            compartment.safe_ambient_pressure(1.01),
            Err(BodyModelError::InvalidGradientFactor)
        );
        assert_eq!(
            compartment.safe_ambient_pressure(f64::NAN),
            Err(BodyModelError::InvalidGradientFactor)
        );
    }
}