From 88079426035862e1299dafa8d06a7ca89583c878 Mon Sep 17 00:00:00 2001
From: Luis Aleixo <luis.aleixo@cern.ch>
Date: Fri, 10 Sep 2021 10:10:44 +0200
Subject: [PATCH] f_dep in function of diameter graph

---
 cara/plot_output.py   | 12 +++++---
 cara/results_paper.py | 66 +++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 74 insertions(+), 4 deletions(-)

diff --git a/cara/plot_output.py b/cara/plot_output.py
index c5dc5fc9..88e51f20 100644
--- a/cara/plot_output.py
+++ b/cara/plot_output.py
@@ -4,7 +4,7 @@ Date: <date>
 Code version: <code version>
 Availability: <where it's located> """
 
-from cara.models import InfectedPopulation
+from cara.models import ExposureModel, InfectedPopulation
 from cara import model_scenarios_paper
 from cara.results_paper import *
 from cara.test_plots import *
@@ -30,12 +30,16 @@ print('\n<<<<<<<<<<< Vlout for Breathing, seated with chosen Cn,B >>>>>>>>>>>')
 print('\n')
 
 ############ Plots with viral loads and emission rates + statistical data ############
-present_vl_er_histograms(activity='Seated', mask='No mask')
-present_vl_er_histograms(activity='Light activity', mask='No mask')
-present_vl_er_histograms(activity='Heavy exercise', mask='No mask')
+#present_vl_er_histograms(activity='Seated', mask='No mask')
+#present_vl_er_histograms(activity='Light activity', mask='No mask')
+#present_vl_er_histograms(activity='Heavy exercise', mask='No mask')
 
 ############ CDFs for comparing the QR-Values in different scenarios ############
 #generate_cdf_curves()
 
+############ Deposition Fraction Graph ############
+print('\n<<<<<<<<<<< Deposition Fraction for Breathing, seated >>>>>>>>>>>') 
+calculate_deposition_factor()
+
 ############ Used for testing ############
 #exposure_model_from_vl_talking_new_points()
diff --git a/cara/results_paper.py b/cara/results_paper.py
index 2ac06fa2..54f3e8fe 100644
--- a/cara/results_paper.py
+++ b/cara/results_paper.py
@@ -11,6 +11,7 @@ import matplotlib.pyplot as plt
 import pandas as pd
 import matplotlib.lines as mlines
 import matplotlib as mpl
+import math
 
 
 ######### Plot material #########
@@ -415,6 +416,71 @@ def generate_cdf_curves():
     fig.set_figwidth(5)
     plt.show()
 
+############ Deposition Fraction Graph ############ 
+# :D particle diameter - varies between 0.3*10**(-6) and 30*10**(-6)
+# :Pp mass density - 1000 kg/m-3
+# :Uair dynamic viscosity of air - 1.8*10**(-5)
+# :FRC functional residual capacity of lungs - 0.003m3
+# :Vt tidal volume - 0.0004m3
+# :g gravity of earth - 9.8m/s**(-2)
+# :BRk mask - inhalation rate
+# Cs(D) Cunningham slip factor
+#   lambda = 0.065
+#   A1 = 1.246
+#   A2 = 0.42
+#   A3 = 0.87
+
+def calculate_cunningham_slip_factor(d:int):
+    _lambda=0.065
+    A1 = 1.246
+    A2 = 0.42
+    A3 = 0.87
+
+    return 1 + (((2*_lambda)/d) * (A1+(A2*(math.e**((-A3*d)/_lambda)))))
+
+# Deposition factor for a breathing model
+def calculate_deposition_factor():
+
+    exposure_mc = breathing_exposure(activity='Seated', mask='No mask')
+    exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
+
+    Pp = 1000
+    Uair = 1.8*10**-5
+    FRC = 0.003
+    Vt = 0.0004
+    g = 9.8
+    BRk = exposure_model.exposed.activity.inhalation_rate
+
+    diameters = np.linspace(0.3, 30, 200) #particle diameter (multiply later by 10**(-6))
+    
+    fractions = []
+    for d in diameters:
+        d = d*10**(-6)
+        cunningham_slip_factor = calculate_cunningham_slip_factor(d)
+        f_dep = 0.08 + 0.92 / (
+            1 + (4.09*10**-6 * (
+                (((cunningham_slip_factor*Pp*d**2*(BRk/3600))/Uair*FRC)**0.8) + (
+                    0.01*(
+                        ((cunningham_slip_factor*g*Pp*d**2*FRC**(2/3))/(Uair*(BRk/3600))**0.4) * (
+                            (Vt/FRC)**0.8
+                        ) 
+                    )
+                )
+            )**(-2.06)
+        ))
+
+        fractions.append(f_dep)
+    
+    fig = plt.figure()
+    ax = fig.add_subplot(1, 1, 1)
+    ax.plot(diameters, fractions)
+    print(diameters)
+
+    plt.ylabel(
+        'Deposition fraction (f$_{dep}$)', fontsize=14)
+    plt.xlabel('Particle diameter (μm)', fontsize=14)
+
+    plt.show()
 
 ######### Auxiliar functions #########