-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathval_Earth_v.py
61 lines (55 loc) · 2.37 KB
/
val_Earth_v.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
################################################################
# calculate and make LoWo21 Figure S4
# validate raindrop velocity against Earth obs and empirical
# relationships
################################################################
import src.drop_prop as drop_prop
import numpy as np
import matplotlib.pyplot as plt
from src.planet import Planet
# set up Earth planet
# following Beard (1976)
X = np.zeros(5) # composition
X[2] = 0.8 # f_N2 [mol/mol]
X[3] = 0.2 # f_O2 [mol/mol]
T_surf = 293.15 # [K]
p_surf = 1.01325e5 # [Pa]
RH_surf = 0.5 # [ ]
R_p = 1. # [R_earth]
M_p = 1. # [M_earth]
pl = Planet(R_p,T_surf,p_surf,X,'h2o',RH_surf,M_p)
pl.z2x4drdz(0.) # set conditions to surface
n = 100 # number of r values
rs = np.logspace(-6,np.log10(7.25e-3/2.),n) # [m]
# set up empty arrays for velocities calculated using different methods
v_B = np.zeros(n) # Beard (1976)
v_Lorenz = np.zeros(n)
v_Lorenz2 = np.zeros(n)
v_Holzer = np.zeros(n)
v_Salman = np.zeros(n)
v_Ganser_Newton = np.zeros(n)
v_sphere = np.zeros(n)
v_Loth = np.zeros(n)
# calculate v(r) using different methods, mainly different C_D parameterizations
for i,r in enumerate(rs):
v_B[i] = drop_prop.calc_v_Earth_Beard(r,pl)
v_Lorenz[i] = drop_prop.calc_v(r,pl,pl.f_rat,drop_prop.calc_C_D_Lorenz)
v_Lorenz2[i] = drop_prop.calc_v(r,pl,drop_prop.calc_rat_Lorenz,drop_prop.calc_C_D_Lorenz)
v_Holzer[i] = drop_prop.calc_v(r,pl,pl.f_rat,drop_prop.calc_C_D_Holzer)
v_Salman[i] = drop_prop.calc_v(r,pl,pl.f_rat,drop_prop.calc_C_D_Salman)
v_Ganser_Newton[i] = drop_prop.calc_v(r,pl,pl.f_rat,drop_prop.calc_C_D_Ganser_Newton)
v_sphere[i] = drop_prop.calc_v(r,pl,pl.f_rat,drop_prop.calc_C_D_sphere)
v_Loth[i] = drop_prop.calc_v(r,pl,pl.f_rat,drop_prop.calc_C_D_Loth)
# load in experimental data
gunn49 = np.genfromtxt('data/gunn1949_r_m_v_mpers.csv',delimiter=',')
best50 = np.genfromtxt('data/best1950.csv',delimiter=',')
# make plot
plt.plot(rs*1e3,v_Loth,c='indigo',lw=3,label='this work')
plt.plot(rs*1e3,v_B,c='plum',ls='--',lw=2,label='Beard (1976)')
plt.scatter(gunn49[:,0]*1e3,gunn49[:,1],c='plum',label='Gunn & Kinzer (1949)',marker='1',zorder=10)
plt.scatter(best50[:,0]/2.,best50[:,1],c='plum',label='Best (1950)',marker='+',zorder=10)
plt.legend()
plt.ylabel(r'$v_\mathrm{T}$ [m s$^{-1}$]')
plt.xlabel(r'$r_\mathrm{eq}$ [mm]')
plt.savefig('sfigs/sfig04.pdf',transparent=True,bbox_inches='tight')
plt.close()