update to documentation on dilution factor

caimira/docs/full_diameter_dependence.rst

@@ -139,8 +139,24 @@ Very similar to what we did with the **emission rate**, we need to calculate the
 During a given exposure time, multiple short-range interactions can be defined in the model.
 In addition, for each individual interaction, the expiration type may be different.
 
-To calculate the short-range component, we first need to calculate what is the **concentration at the jet origin**, that depends on the diameter :math:`D`. 
-The initial concentration of virions at the mouth/nose, :math:`C_{0, \mathrm{SR}}(D)` is calculated as follows:
+To calculate the short-range component, we first need to calculate what is the **dilution factor**, that depends on the distance :math:`x` as a random variable, from a log normal distribution in :meth:`caimira.monte_carlo.data.short_range_distances`.
+This factor is calculated in a two-stage model, based in a transition point, defined as follows:
+
+:math:`\mathrm{xstar}=𝛽_{\mathrm{x1}} \cdot (\mathrm{BR} \cdot u_{0})^\frac{1}{4} \cdot (\mathrm{tstar} + t_{0})^\frac{1}{2} - x_{0}`,
+
+where :math:`\mathrm{BR}` is the expired flow rate during the expiration period, converted from :math:`m^{3} h^{-1}` to :math:`m^{3} s^{-1}`, :math:`u_{0}` is the expired jet speed (in :math:`m s^{-1}`) given by :math:`u_{0}=\frac{\mathrm{BR}}{A_{m}}`, being :math:`A_{m}` the area of the mouth assuming a perfect circle (average mouth diameter :math:`D` of `0.02m`).
+The time of the transition point :math:`\mathrm{tstar}` is defined as 2s. The distance of virtual origin :math:`x_{0}` given by :math:`x_{0}=\frac{D}{2𝛽_{\mathrm{r1}}}` (in m), and the time of virtual origin on puff-like stage is given by :math:`t_{0}=(\frac{x_{0}}{𝛽_{\mathrm{x1}}})^2 \cdot (\mathrm{BR} \cdot u_{0})^-\frac{1}{2}` (in s).
+Having the distance for the transition point, we can calculate the dilution factor at the transition point, defined as follows:
+
+:math:`\mathrm{Sxstar}=2𝛽_{\mathrm{r1}}\frac{(xstar + x_{0})}{D}`.
+
+The remaining dilution factors, either in the jet- or puff-like stages are calculated as follows:
+
+:math:`\mathrm{factors}(x)=\begin{cases}\hfil 2𝛽_{\mathrm{r1}}\frac{(x + x_{0})}{D} & \textrm{if } x < \mathrm{xstar},\\\hfil \mathrm{Sxstar} \cdot \biggl(1 + \frac{𝛽_{\mathrm{r2}}(x - xstar)}{𝛽_{\mathrm{r1}}(xstar + x_{0})}\biggl)^3 & \textrm{if } x > \mathrm{xstar}.\end{cases}`
+
+The variables :math:`𝛽_{\mathrm{r1}}`, :math:`𝛽_{\mathrm{r2}}` and :math:`𝛽_{\mathrm{x1}}` are defined as `0.18`, `0.2`, and `2.4` respectively. The dilution factor for each distance `x` is then stored in the :math:`\mathrm{factors}` array that is returned by the method.
+
+Having the dilution factors, the **initial concentration of virions at the mouth/nose**, :math:`C_{0, \mathrm{SR}}(D)`, is calculated as follows:
 
 :math:`C_{0, \mathrm{SR}}(D) = N_p(D) \cdot V_p(D) \cdot \mathrm{vl_{in}} \cdot 10^{-6}`, 
 given by :meth:`caimira.models.Expiration.jet_origin_concentration`. It computes the same quantity as :meth:`caimira.models.Expiration.aerosols`, except for the mask inclusion. As previously mentioned, it is normalized by the **viral load**, which is a diameter-independent property. 
@@ -155,7 +171,7 @@ The former operation is given in method :meth:`caimira.models.ShortRangeModel._l
 one solution would be to recompute the values a second time using :math:`D_{\mathrm{max}} = 100\mathrm{μm}`;
 or perform a approximation using linear interpolation, which is possible and more effective in terms of performance. We decided to adopt the interpolation solution.
 The set of points with a known value are given by the default expiration particle diameters for long-range, i.e. from 0 to 30 :math:`\mathrm{μm}`.
-The set of points we want the interpolated values are given by the short-range expiration particle diameters, i.e. from 0 to 100:math:`\mathrm{μm}`. 
+The set of points we want the interpolated values are given by the short-range expiration particle diameters, i.e. from 0 to 100 :math:`\mathrm{μm}`. 
 
 To summarize, in the code, :math:`C_{\mathrm{SR}}(t, D)` is computed as follows:
 

caimira/models.py

@@ -1156,7 +1156,7 @@ class ShortRangeModel:
         '''
         # Average mouth diameter
         D = 0.02
-        # Convert Breathing rate from m3/h to m3/s
+        # The expired flow rate during the expiration period, converted from m3/h to m3/s
         BR = np.array(self.activity.exhalation_rate/3600.)
         # Area of the mouth assuming a perfect circle
         Am = np.pi*(D**2)/4 
@@ -1164,30 +1164,28 @@ class ShortRangeModel:
         u0 = np.array(BR/Am)
 
         tstar = 2.0
-        Cr1 = 0.18
-        Cr2 = 0.2
-        Cx1 = 2.4
+        𝛽r1 = 0.18
+        𝛽r2 = 0.2
+        𝛽x1 = 2.4
 
-        # The expired flow rate during the expiration period, m^3/s
-        Q0 = u0 * np.pi/4*D**2 
         # Parameters in the jet-like stage
-        x01 = D/2/Cr1
+        x0 = D/2/𝛽r1
         # Time of virtual origin
-        t01 = (x01/Cx1)**2 * (Q0*u0)**(-0.5)
+        t0 = (x0/𝛽x1)**2 * (BR*u0)**(-0.5)
         # The transition point, m
-        xstar = np.array(Cx1*(Q0*u0)**0.25*(tstar + t01)**0.5 - x01)
+        xstar = np.array(𝛽x1*(BR*u0)**0.25*(tstar + t0)**0.5 - x0)
         # Dilution factor at the transition point xstar
-        Sxstar = np.array(2*Cr1*(xstar+x01)/D)
+        Sxstar = np.array(2*𝛽r1*(xstar+x0)/D)
 
         distances = np.array(self.distance)
 
         factors = np.empty(distances.shape, dtype=np.float64)
-        factors[distances < xstar] = 2*Cr1*(distances[distances < xstar]
-                                        + x01)/D
+        factors[distances < xstar] = 2*𝛽r1*(distances[distances < xstar]
+                                        + x0)/D
         factors[distances >= xstar] = Sxstar[distances >= xstar]*(1 +
-            Cr2*(distances[distances >= xstar] -
-            xstar[distances >= xstar])/Cr1/(xstar[distances >= xstar]
-            + x01))**3
+            𝛽r2*(distances[distances >= xstar] -
+            xstar[distances >= xstar])/𝛽r1/(xstar[distances >= xstar]
+            + x0))**3
         return factors
 
     def _long_range_normed_concentration(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat: