-
Notifications
You must be signed in to change notification settings - Fork 1
/
val_Earth_shape.py
65 lines (55 loc) · 2.65 KB
/
val_Earth_shape.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
62
63
64
65
################################################################
# calculate and make LoWo21 Figure S3
# validate raindrop shape ratio against Earth obs and models
################################################################
import numpy as np
import matplotlib.pyplot as plt
import src.drop_prop as drop_prop
from src.planet import Planet
# set up Earth planet
X = np.zeros(5) # composition
X[2] = 0.79 # f_N2 [mol/mol]
X[3] = 0.21 # f_O2 [mol/mol]
T_surf = 290. # [K]
p_surf = 1.01325e5 # [Pa]
RH_surf = 0.75 # [ ]
R_p = 1. # [R_earth]
M_p = 1. # [M_earth]
Earth = Planet(R_p,T_surf,p_surf,X,'h2o',RH_surf,M_p)
# maximum stable raindrop size
# from Rayleigh-Taylor instability, ℓ_RT = 0.5πa
max_r_Earth = drop_prop.calc_r_max_RT(Earth,np.pi/2.)
n = 100 # number of points
rs = np.linspace(0.01e-3,max_r_Earth,n) # [m] rs to calculate ratio for
rat = np.zeros(n) # [ ] array to fill with b/a shape ratio
for i,r in enumerate(rs):
rat[i] = Earth.f_rat(r,Earth) # calc shape ratio using Beard method
# ratios vs r from other sources
# load in rat vs r data
thurai_2dvd = np.genfromtxt('data/Thurai2009_2DVD_a_v_d_mm.dat')
thurai_tunnel = np.genfromtxt('data/Thurai2009_windtunnel_a_v_d_mm.dat')
beard1991 = np.genfromtxt('data/Beard1991_a_v_d_mm.dat')
# Pruppacher & Beard (1970)
r_lim = np.linspace(0.5,max_r_Earth*1e3,n) # [mm]
rat_PB70 = 1.030 - 0.062*2*r_lim # [ ]
# Beard & Chuang (1987)
r_BC87 = np.array([1.,2.,3.,4,5.,6.,7.,8.])/2. # [mm]
rat_BC87 = np.array([0.983,0.928,0.853,0.778,0.708,0.642,0.581,0.521]) # [ ]
# Pruppacher & Pitt (1971)
r_PP71 = np.array([0.0170,0.0305,0.0433,0.0532,0.0620,0.11,0.14,0.15,0.18,0.20,0.25,0.29,0.30,0.35,0.40])*10. # [mm]
rat_PP71 = np.array([0.9993,0.9959,0.9892,0.9813,0.9735,0.916,0.865,0.847,0.795,0.762,0.701,0.664,0.655,0.621,0.583]) # [ ]
# make plot
plt.plot(rs*1e3, rat, c='indigo', lw=2,label='our model',zorder=8) # convert rs from m to mm
plt.plot(r_lim,rat_PB70,c='0.75',lw=2,ls='--',label='Pruppacher & Beard (1970)')
plt.scatter(r_PP71,rat_PP71,marker='o',label='Pruppacher & Pitt (1971)',facecolors='none',edgecolor='0.25',zorder=9)
plt.scatter(r_BC87,rat_BC87,facecolor='0.25',marker='x',label='Beard & Chuang (1987)',zorder=10)
plt.scatter(beard1991[:,0]/2,beard1991[:,1],c='plum',zorder=9,marker='_', label='Beard et al. (1991)')
plt.scatter(thurai_2dvd[:-7,0]/2,thurai_2dvd[:-7,1],c='plum',marker='1',zorder=9,label='Thurai et al. (2009)')
# :-7 cuts off points beyond r_max
plt.scatter(thurai_tunnel[:,0]/2,thurai_tunnel[:,1],c='plum',zorder=9,marker='1')
plt.xlabel(r'$r_\mathrm{eq}$ [mm]')
plt.ylabel(r'$b/a$ [ ]')
plt.ylim(0,1.05)
plt.legend()
plt.savefig('sfigs/sfig03.pdf',bbox_inches='tight',transparent=True)
plt.close()