From e7a1c831b723209371e937a22eee44c9760c4ec1 Mon Sep 17 00:00:00 2001
From: Luis Aleixo <luis.aleixo@cern.ch>
Date: Tue, 20 Jul 2021 12:36:21 +0200
Subject: [PATCH] Concentration profile plot with cumulative dose #fixes164

---
 cara/apps/calculator/report_generator.py | 38 ++++++++++++++++++------
 1 file changed, 29 insertions(+), 9 deletions(-)

diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py
index a899a18f..4d2b8dc0 100644
--- a/cara/apps/calculator/report_generator.py
+++ b/cara/apps/calculator/report_generator.py
@@ -105,18 +105,20 @@ def img2base64(img_data) -> str:
     # A src suitable for a tag such as f'<img id="scenario_concentration_plot" src="{result}">.
     return f'data:image/png;base64,{pic_hash}'
 
-
 def plot(times, concentrations, model: models.ExposureModel):
     fig = plt.figure()
     ax = fig.add_subplot(1, 1, 1)
     datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) for time in times]
-    ax.plot(datetimes, concentrations, lw=2, color='#1f77b4', label='Mean concentration')
+    
+    #Concentration as mean viral concentration (virion m$^{-3}$)
+    concentrations = [c * model.concentration_model.virus.quantum_infectious_dose for c in concentrations]
+    
+    ax.plot(datetimes, concentrations, lw=2, color='#1f77b4', label='Mean viral concentration')
     ax.spines['right'].set_visible(False)
-    ax.spines['top'].set_visible(False)
 
-    ax.set_xlabel('Time of day')
-    ax.set_ylabel('Mean concentration ($q/m^3$)')
-    ax.set_title('Mean concentration of infectious quanta')
+    ax.set_xlabel('Time of day', fontsize=14)
+    ax.set_ylabel('Mean viral concentration\n(virion m$^{-3}$)', fontsize=14)
+    ax.set_title('Concentration profile')
     ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M"))
 
     # Plot presence of exposed person
@@ -128,13 +130,31 @@ def plot(times, concentrations, model: models.ExposureModel):
             label="Presence of exposed person(s)" if i == 0 else ""
         )
 
+    #See CERN-OPEN-2021-004, p. 15, eq. 16. - Cumulative Dose
+    factor = 0.6 * np.mean(model.exposed.activity.inhalation_rate) * (1 - model.exposed.mask.η_inhale)
+    present_indexes = np.array([model.exposed.person_present(t) for t in times])
+    modified_concentrations = np.array(concentrations)
+    modified_concentrations[~present_indexes] = 0
+    
+    qds = [np.trapz(modified_concentrations[:i + 1], times[:i + 1]) * factor for i in range(len(times))]
+    
+    ax1 = ax.twinx()
+    ax1.plot(datetimes, qds, label='qD - Mean viral concentration', color='#1f77b4', linestyle='dotted')
+    ax1.spines["right"].set_linestyle("--")
+    ax1.spines["right"].set_linestyle((0,(1,5)))
+    ax1.set_ylabel('Mean cumulative dose\n(virion)', fontsize=14)
+    ax1.set_ylim(ax1.get_ylim()[0], ax1.get_ylim()[1])
+    ax1.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M"))
+
     # Place a legend outside of the axes itself.
-    ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
-    ax.set_ylim(0)
+    fig.legend(bbox_to_anchor=(1.05, 0.9), loc='upper left')
+    ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1])
+    # Remove top spines
+    ax.spines['top'].set_visible(False)
+    ax1.spines['top'].set_visible(False)
 
     return fig
 
-
 def minutes_to_time(minutes: int) -> str:
     minute_string = str(minutes % 60)
     minute_string = "0" * (2 - len(minute_string)) + minute_string