From 0ce14e56f3f33a1f68cebe41be1f1f2525c03fbb Mon Sep 17 00:00:00 2001
From: Andrejh <andre.henriques@cern.ch>
Date: Sun, 12 Sep 2021 14:27:49 +0200
Subject: [PATCH] add results for shouting

---
 cara/model_scenarios_paper.py | 31 +++++++++++++++++++++++++++++++
 cara/plot_output.py           |  6 +++++-
 cara/results_paper.py         | 33 +++++++++++++++++++++++++++++++++
 3 files changed, 69 insertions(+), 1 deletion(-)

diff --git a/cara/model_scenarios_paper.py b/cara/model_scenarios_paper.py
index 901db116..608d86ea 100644
--- a/cara/model_scenarios_paper.py
+++ b/cara/model_scenarios_paper.py
@@ -183,6 +183,37 @@ def talking_exposure_vl(vl):
     )
     return exposure_mc
 
+######### Shouting model for specific viral load ###########
+def shouting_exposure_vl(vl):
+    exposure_mc = mc.ExposureModel(
+        concentration_model=mc.ConcentrationModel(
+            room=models.Room(volume=100, humidity=0.5),
+            ventilation=models.AirChange(
+                active=models.SpecificInterval(((0, 24),)),
+                air_exch=0.25,
+            ),
+            infected=mc.InfectedPopulation(
+                number=1,
+                virus=mc.Virus(
+                    viral_load_in_sputum=10**vl,
+                    infectious_dose=50.,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types['No mask'],
+                activity=activity_distributions['Light activity'],
+                expiration=models.Expiration.types['Shouting'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types['Light activity'],
+            mask=models.Mask.types["No mask"],
+        ),
+    )
+    return exposure_mc
+
 ######### Used for CDF Models ###########
 ######### Breathing Models #########
 def breathing_seated_exposure():
diff --git a/cara/plot_output.py b/cara/plot_output.py
index 88e51f20..58e39b38 100644
--- a/cara/plot_output.py
+++ b/cara/plot_output.py
@@ -18,7 +18,11 @@ print('\n<<<<<<<<<<< Vlout for Talking, seated >>>>>>>>>>>')
 
 # Exhaled virions while breathing, seated #
 print('\n<<<<<<<<<<< Vlout for Breathing, seated >>>>>>>>>>>')
-#exposure_model_from_vl_breathing()
+exposure_model_from_vl_breathing()
+
+# Exhaled virions while breathing, light activity #
+print('\n<<<<<<<<<<< Vlout for Shouting, light activity >>>>>>>>>>>')
+exposure_model_from_vl_shouting()
 
 # Exhaled virions while talking according to BLO model, seated #
 print('\n<<<<<<<<<<< Vlout for Talking, seated with chosen Cn,L >>>>>>>>>>>')
diff --git a/cara/results_paper.py b/cara/results_paper.py
index a94f4df5..75cdabbc 100644
--- a/cara/results_paper.py
+++ b/cara/results_paper.py
@@ -274,6 +274,39 @@ def exposure_model_from_vl_talking_cn():
     plt.xlabel('NP viral load, $\mathrm{vl_{in}}$\n(RNA copies)', fontsize=14)
     plt.show()
 
+####### Shouting ########
+
+
+def exposure_model_from_vl_shouting():
+    fig = plt.figure()
+    ax = fig.add_subplot(1, 1, 1)
+
+    er_means_1h = []
+
+    for vl in tqdm(viral_loads):
+        exposure_mc = shouting_exposure_vl(vl)
+        exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
+        emission_rate_1h = exposure_model.concentration_model.infected.emission_rate_when_present(cn_B=0.06, cn_L=0.2)
+        er_means_1h.append(np.mean(emission_rate_1h))
+
+    ax.plot(viral_loads, er_means_1h)
+    ax.set_yscale('log')
+
+    ratio_1h = np.mean(10**viral_loads / er_means_1h)
+    print('Mean swab-to-aersol vl ratio emission rate per hour:')
+    print(format(ratio_1h, "5.1e"))
+
+    ############ Plot ############
+    plt.title('',
+              fontsize=16, fontweight="bold")
+    plt.ylabel(
+        'Aerosol viral load, $\mathrm{vl_{out}}$\n(RNA copies)', fontsize=14)
+    plt.xticks(ticks=[i for i in range(2, 13)], labels=[
+        '$10^{' + str(i) + '}$' for i in range(2, 13)])
+    plt.xlabel('NP viral load, $\mathrm{vl_{in}}$\n(RNA copies)', fontsize=14)
+    plt.show()
+
+
 ############ Plots with viral loads and emission rates ############
 ############ Statistical Data ############