diff --git a/cara/models.py b/cara/models.py index a55ccf97..71c532d3 100644 --- a/cara/models.py +++ b/cara/models.py @@ -552,56 +552,20 @@ class _ExpirationBase: @dataclass(frozen=True) class Expiration(_ExpirationBase): """ - BLO model for the expiration (G. Johnson et al., Modality of human - expired aerosol size distributions, Journal of Aerosol Science, - vol. 42, no. 12, pp. 839 – 851, 2011, - https://doi.org/10.1016/j.jaerosci.2011.07.009). - Here all diameters (d) are in microns. + Model for the expiration. For a given diameter of aerosol, provides + the aerosol volume, weighted by the mask outward efficiency when + applicable. """ - #: factors assigned to resp. the B, L and O modes. They are - # charateristics of the kind of expiratory activity (e.g. breathing, - # speaking, singing, or shouting). - BLO_factors: typing.Tuple[float, float, float] - #: diameter of the aerosol in microns diameter: _VectorisedFloat - @cached() - def aerosols_(self, mask: Mask): - """ Result is in mL.cm^-3 """ - def volume(d): - return (np.pi * d**3) / 6. - - def _Bmode(d: float) -> float: - # B-mode (see ref. above). - return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.262364)) * - np.exp(-1 * (np.log(d) - 0.989541) ** 2 / (2 * 0.262364 ** 2))) - - def _Lmode(d: float) -> float: - # L-mode (see ref. above). - return ( (1 / d) * (1.0 / (np.sqrt(2 * np.pi) * 0.506818)) * - np.exp(-1 * (np.log(d) - 1.38629) ** 2 / (2 * 0.506818 ** 2))) - - def _Omode(d: float) -> float: - # O-mode (see ref. above). - return ( (1 / d) * (0.0010008 / (np.sqrt(2 * np.pi) * 0.585005)) * - np.exp(-1 * (np.log(d) - 4.97673) ** 2 / (2 * 0.585005 ** 2))) - - def integrand(d: float) -> float: - return (self.BLO_factors[0] * _Bmode(d) + - self.BLO_factors[1] * _Lmode(d) + - self.BLO_factors[2] * _Omode(d) - ) * volume(d) * (1 - mask.exhale_efficiency(d)) - - # final result converted from microns^3/cm3 to mL/cm^3 - return scipy.integrate.quad(integrand, 0.1, 30.)[0]*1e-12 - @cached() def aerosols(self, mask: Mask): """ Result is in mL.cm^-3 """ def volume(d): return (np.pi * d**3) / 6. + # final result converted from microns^3/cm3 to mL/cm^3 return (volume(self.diameter) * (1 - mask.exhale_efficiency(self.diameter))) * 1e-12 @@ -632,12 +596,14 @@ class MultipleExpiration(_ExpirationBase): # Typical expirations. The aerosol diameter given is an equivalent # diameter, chosen in such a way that the aerosol volume is -# the same as the total aerosol volume given by the full BLO model. +# the same as the total aerosol volume given by the full BLO model +# (integrated between 0.1 and 30 microns) +# The correspondence with the BlO coefficients is given. _ExpirationBase.types = { - 'Breathing': Expiration((1., 0., 0.)), - 'Talking': Expiration((1., 1., 1.)), - 'Shouting': Expiration((1., 5., 5.)), - 'Singing': Expiration((1., 5., 5.)), + 'Breathing': Expiration(1.3844), # corresponds to B/L/O coefficients of (1, 0, 0) + 'Talking': Expiration(5.8925), # corresponds to B/L/O coefficients of (1, 1, 1) + 'Shouting': Expiration(10.0411), # corresponds to B/L/O coefficients of (1, 5, 5) + 'Singing': Expiration(10.0411), # corresponds to B/L/O coefficients of (1, 5, 5) } diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py index f82bbd63..e3497efe 100644 --- a/cara/tests/test_expiration.py +++ b/cara/tests/test_expiration.py @@ -9,7 +9,7 @@ from cara import models def test_multiple_wrong_weight_size(): weights = (1., 2., 3.) - e_base = models.Expiration((1., 0., 0.), 2.5) + e_base = models.Expiration(2.5) with pytest.raises( ValueError, match=re.escape("expirations and weigths should contain the" @@ -20,12 +20,12 @@ def test_multiple_wrong_weight_size(): def test_multiple(): weights = (1., 1.) - e1 = models.Expiration((1., 0., 0.), 2.5) - e2 = models.Expiration((4., 5., 5.), 2.5) - e_expected = models.Expiration((2.5, 2.5, 2.5), 2.5) - e = models.MultipleExpiration([e1, e2], weights) mask = models.Mask.types['Type I'] - npt.assert_almost_equal(e_expected.aerosols(mask), e.aerosols(mask)) + e1 = models.Expiration.types['Breathing'] + e2 = models.Expiration.types['Singing'] + aerosol_expected = (e1.aerosols(mask) + e2.aerosols(mask))/2. + e = models.MultipleExpiration([e1, e2], weights) + npt.assert_almost_equal(aerosol_expected, e.aerosols(mask)) # expected values obtained from another code diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index f153d389..fbf9f105 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -218,10 +218,11 @@ def skagit_chorale_mc(): presence=mc.SpecificInterval(((0., 2.5),)), mask=models.Mask.types["No mask"], activity=activity_distributions['Light activity'], - expiration=models.Expiration((5., 5., 5.), 10.0761), + expiration=models.Expiration(10.0761), # The aerosol diameter given (10.0761 microns) is an equivalent # diameter, chosen in such a way that the aerosol volume is - # the same as the total aerosol volume given by the full BLO model. + # the same as the total aerosol volume given by the full BLO model + # with (5, 5, 5) for the B/L/O weights. ), ) return mc.ExposureModel(