-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmie_droplets.py
96 lines (80 loc) · 2.64 KB
/
mie_droplets.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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
"""
Plot the scattering cross section as a function of wavelength for 1 micron water droplets
"""
import numpy as np
import matplotlib.pyplot as plt
import miepython
import urllib
import segelstein
def mie_scatter (radius, lam,m):
x = 2*np.pi*radius/lam
num = len(lam)
qqsca = np.zeros(num)
qqext = np.zeros(num)
for i in range(num) :
qext, qsca, qback, g = miepython.mie(m[i],x[i])
qqsca[i]=qsca
qqext[i]=qext
ssa = qqsca/qqext #single scattering albedo
return qqsca, ssa
#load refractive index of water data
url = 'http://www.philiplaven.com/Segelstein.txt'
lam,rfr,rfi = segelstein.segel(url)
lam = np.array(lam) #microns
rfr = np.array(rfr)
rfi = np.array(rfi)
m = rfr + 1j*rfi
radius1 = 7.5 # in microns
radius2 = 2.5 # in microns
lam1 = lam[np.where(np.logical_and(lam>=0.4, lam <=1.7))] # also in microns
rfr1 = rfr[np.where(np.logical_and(lam >= 0.4, lam <= 1.7))]
rfi1 = rfi[np.where(np.logical_and(lam >= 0.4, lam <= 1.7))]
m1 = rfr1 + 1j*rfi1
lam2 = lam[np.where(np.logical_and(lam>=10, lam <=30))] # also in microns
rfr2 = rfr[np.where(np.logical_and(lam >= 10, lam <= 30))]
rfi2 = rfi[np.where(np.logical_and(lam >= 10, lam <= 30))]
m2 = rfr2 + 1j*rfi2
qe_a1, ssa_a1 = mie_scatter(radius1, lam1,m1)
qe_a2, ssa_a2 = mie_scatter(radius2, lam1,m1)
qe_b1, ssa_b1 = mie_scatter(radius1, lam2,m2)
qe_b2, ssa_b2 = mie_scatter(radius2, lam2,m2)
plt.figure(figsize = (8,6))
plt.plot(lam1*1000,qe_a1,'k--',label='d = 5 micron')
plt.plot(lam1*1000,qe_a2,'r--',label='d = 15 micron')
plt.title(r"Mie scattering efficiency - 2a")
plt.xlabel("Wavelength (nm)")
plt.ylabel("Scattering Efficiency")
#plt.ylim(1,3)
plt.legend()
plt.savefig('2.jpg')
plt.show()
plt.figure(figsize = (8,6))
plt.plot(lam1*1000,ssa_a1,'k--',label='d = 5 micron')
plt.plot(lam1*1000,ssa_a2,'r--',label='d = 15 micron')
plt.title(r"Single scatterig albedo- 2a")
plt.xlabel("Wavelength (nm)")
plt.ylabel("Single Scattering Albedo")
plt.legend()
plt.savefig('2b.jpg')
#plt.ylim(1,3)
plt.show()
plt.figure(figsize = (8,6))
plt.plot(lam2*1000,qe_b1,'k--',label='d = 5 micron')
plt.plot(lam2*1000,qe_b2,'r--',label='d = 15 micron')
plt.title(r"Mie scattering efficiency - 2b")
plt.xlabel("Wavelength (nm)")
plt.ylabel("Scattering Efficiency")
plt.legend()
plt.savefig('2c.jpg')
#plt.ylim(1,3)
plt.show()
plt.figure(figsize = (8,6))
plt.plot(lam2*1000,ssa_b1,'k--',label='d = 5 micron')
plt.plot(lam2*1000,ssa_b2,'r--',label='d = 15 micron')
plt.title(r"Single scattering albedo- 2b")
plt.xlabel("Wavelength (nm)")
plt.ylabel("Single Scattering Albedo")
plt.legend()
plt.savefig('2d.jpg')
#plt.ylim(1,3)
plt.show()