From 9810ee3f1fa1dacf4b0ff0372638154eca6e021a Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 15 Sep 2021 10:38:17 +0200 Subject: [PATCH] deposition_fraction_plot with correct values --- cara/results_paper.py | 123 +++++++++++++----------------------------- 1 file changed, 38 insertions(+), 85 deletions(-) diff --git a/cara/results_paper.py b/cara/results_paper.py index d85d5234..0658940c 100644 --- a/cara/results_paper.py +++ b/cara/results_paper.py @@ -470,107 +470,60 @@ def generate_cdf_curves(): 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*10**-6 - A1 = 1.246 - A2 = 0.42 - A3 = 0.87 - - return 1 + (((2*_lambda)/d) * (A1+(A2*(math.e**((-A3*d)/_lambda))))) - -# Deposition factors for a breathing model +# def calculate_deposition_factor(): fig = plt.figure() ax = fig.add_subplot(1, 1, 1) - - ############ Breathing models ############ - br_seated = breathing_seated_exposure() - br_seated_model = br_seated.build_model(size=SAMPLE_SIZE) - - br_light_activity = breathing_light_activity_exposure() - br_light_activity_model = br_light_activity.build_model(size=SAMPLE_SIZE) - - br_heavy_exercise = breathing_heavy_exercise_exposure() - br_heavy_exercise_model = br_heavy_exercise.build_model(size=SAMPLE_SIZE) - - rho_p = 1000 - mu_air = 1.8*10**-5 - FRC = 0.003 - Vt = 0.0004 - g = 9.81 - k = 1.38*10**-23 - T = 300 - diameters = np.linspace(0.001, 0.01, 1000) #particle diameter (multiply later by 10**(-6)) - diameters = np.append(diameters, np.linspace(0.01, 100, 400)) - activity_fractions = [] - for scenario in (br_seated_model, br_light_activity_model, br_heavy_exercise_model): - BRk = scenario.exposed.activity.inhalation_rate - fractions = [] - for d_μm in diameters: - d = d_μm*10**(-6) - cunningham_slip_factor = calculate_cunningham_slip_factor(d) + diameters = np.linspace(0.001, 0.01, 1000) #particle diameter (μm) + diameters = np.append(diameters, np.linspace(0.01, 0.1, 1000)) + diameters = np.append(diameters, np.linspace(0.1, 1., 1000)) + diameters = np.append(diameters, np.linspace(1., 10., 1000)) + diameters = np.append(diameters, np.linspace(10, 100, 1000)) - stk_ = (cunningham_slip_factor*rho_p*(d**2)*(BRk/3600)) / (mu_air * FRC) - t_ = (cunningham_slip_factor*g*rho_p*(d**2)*(FRC**(2/3))) / (mu_air * (BRk / 3600)) + fractions_et = [] + fractions_tb = [] + fractions_al = [] + fractions_df = [] + for d in diameters: + + IF = 1 - 0.5 * (1 - (1 / (1 + (0.00076*(d**2.8))))) + DF_et = IF * ( + (1 / ( 1 + np.exp(6.84 + 1.183 * np.log(d)))) + (1 / (1 + np.exp(0.924 - 1.885 * np.log(d)))) + ) + fractions_et.append(DF_et) - f_dep_ine = 0.08 + (0.92 / ( - 1 + ((4.09*10**(-6) * ( - (stk_**0.8) + 0.01 * (t_**0.4) * ((Vt/FRC)**0.8) - ) ** -2.06)) - )) + DF_tb = (0.00352/d) * (np.exp(-0.234*((np.log(d) + 3.40)**2)) + (63.9 * np.exp(-0.819*((np.log(d) - 1.61)**2)))) + fractions_tb.append(DF_tb) - X_ = ((k * T * cunningham_slip_factor) / ((3 * math.pi * mu_air * d)) * (((Vt**(1/3)) / (BRk / 3600)))) - y_ = Vt / FRC - - f_dep_diff = 1 - ( - 1 / ((7380 * (X_ **0.539) * y_**0.884) + 1) - ) - f_dep_sed = (0.431 * f_dep_ine) + (0.541 * f_dep_diff) + (1.060 * (f_dep_ine**2)) + (0.685 * (f_dep_diff**2)) - (1.521 * f_dep_ine * f_dep_diff) - - # Eq. S.1 - if d_μm < 0.1: - f_dep = f_dep_diff - elif d_μm >= 0.1 and d_μm < 1.: - f_dep = f_dep_sed - elif d_μm > 1: - f_dep = f_dep_ine + DF_al = (0.0155/d) * (np.exp(-0.416*((np.log(d) + 2.84)**2)) + (19.11*np.exp(-0.482 * ((np.log(d) - 1.362)**2)))) + fractions_al.append(DF_al) - fractions.append(f_dep) - activity_fractions.append(fractions) + DF = IF * (0.0587 + (0.911/(1 + np.exp(4.77 + 1.485 * np.log(d)))) + (0.943/(1 + np.exp(0.508 - 2.58 * np.log(d))))) + fractions_df.append(DF) - # Colors can be changed here - colors = ['skyblue', 'royalblue', 'navy'] - - breathing_rates = ['Seated', 'Light activity', 'Heavy activity'] - - lines = [mlines.Line2D([], [], color=color, markersize=15, label=label) - for color, label in zip(colors, breathing_rates)] - - for i in range(3): - ax.plot(diameters, activity_fractions[i], c=colors[i]) + ax.plot(diameters, fractions_df, label = 'Total') + ax.plot(diameters, fractions_et, label = 'Extrathoracic (ET)', ls='--', lw=0.9) + ax.plot(diameters, fractions_tb, label = 'Tracheobronchial (TB)', ls='--', lw=0.9) + ax.plot(diameters, fractions_al, label = 'Alveolar-interstitial (AI)', ls='--', lw=0.9) ax.grid(linestyle='--') - ax.legend(handles=lines, loc='upper left') ax.set_xscale('log') + ax.margins(x=0, y=0) + plt.legend(bbox_to_anchor=(1.04,1), loc="upper left") + y_ticks = [0., 0.2, 0.4, 0.6, 0.8, 1] + x_ticks = [0.001, 0.01, 0.1, 1, 10, 100] + plt.yticks(ticks=y_ticks, labels=[ + str(i) for i in y_ticks]) + plt.xticks(ticks=x_ticks, labels=[ + str(i) for i in x_ticks]) plt.xlabel('Particle diameter (μm)', fontsize=14) plt.ylabel('Deposition fraction (f$_{dep}$)', fontsize=14) - + + fig.set_figwidth(10) + plt.tight_layout() plt.show() ######### Auxiliar functions #########