Python Mie code
PyMieScatt is a python mie code based on based on Bohren and Huffman's Mie Theory
derivations which calculates extinction efficiency (Qext), scattering efficiency (Qsca),
backscattering efficiency (Qback) and asymmetry parameter (g). It requires the refractive
index of the scattering particle and the size parameter.
PyMieScatt is available at https://github.com/bsumlin/PyMieScatt and can be installed by
running the following command in terminal.
pip install PyMieScatt or conda install -c conda-forge pymiescatt
Requirements
Python environment
PyMieScatt
NumPy
SciPy
matplotlib-pyplot
In [12]: import math # Importing math module for calculations
import PyMieScatt as mie # Importing miepython module for mie calcula
import numpy as np # Importing numpy module for array calculati
from scipy.integrate import trapezoid as trap # Importing scipy module for
# integration
import matplotlib.pyplot as plt # Importing pyplot module for plotting
Single particle
Mie code calculates the mie efficiencies and asymmetry parameter from the particle
refractive index, particle diameter and wavelength of interacting radiation
In [13]: wavelength = 0.5 # Wavelength of incident radiation
radius = 0.1 # Radius of the scattering particle
ref_index = 1.45+0.45j # Refractive index of the scattering particl
# the format n+ik
mie_params = mie.MieQ(ref_index,wavelength,radius*2,asDict=True)
mie_params
{'Qext': 1.625363260122038,
Out[13]:
'Qsca': 0.45142608901140946,
'Qabs': 1.1739371711106286,
'g': 0.35644716992630443,
'Qpr': 1.4644537082630211,
'Qback': 0.19816623330368757,
'Qratio': 0.43897824722017575}
Polydisperse particles
In atmosphere particles exist in different size and composition. In order to assess the
scattering properties of atmospheric particles, we have to take account of their size
distribution. Size distribution of atmospheric particles can be best assumed as lognormal
which can be described by its
mod radius
standard deviation
upper and lower bounds of particle size
2 2
dn(r)/dr = (n/(√2π ln σr) exp [− ln (r/rm )/2(lnσ) ]
In [52]: n = 100 # Number of particles
rm = 0.0118 # Mod radius
sigma = 2 # Standard deviation
rmin = 0.005 # Minimum radius
rmax = 10 # Maximum radius
rad_array = np.arange(rmin,rmax,0.001)
dnr_array = []
for r in rad_array:
dnr = math.sqrt(math.pi/2.0)*(n/np.log(sigma))*r*math.exp(-0.5*(((np.log
dnr_array.append(dnr)
fig = plt.Figure(figsize=(5,5))
plt.plot(rad_array,dnr_array)
plt.xscale('log')
plt.xlabel('radius (um)')
plt.ylabel('dn(r)/dr')
plt.title('Number-size lognormal distribution')
Text(0.5, 1.0, 'Number-size lognormal distribution')
Out[52]:
In [32]: wavelength = 0.5 # Wavelength of incident radiation
ref_index = 1.45+0.45j # Refractive index of the scattering particl
# the format n+ik
n = 1 # Number of particles
rm = 0.0118 # Mod radius
sigma = 2 # Standard deviation
rmin = 0.005 # Minimum radius
rmax = 10 # Maximum radius
rad_array = np.arange(rmin,rmax,0.001)
def miecoeff(n,sigma,rad_array,lamb,rm,m):
bext_array = [] # Arrays to store f(x)
bsca_array = [] # of each radius value
babs_array = [] # for each radius value
g_array = []
for r in rad_array:
mie_params = mie.MieQ(ref_index,wavelength,radius*2,asDict=True) #
# mie parameters
bext = math.sqrt(math.pi/2.0)*(n/np.log(sigma))*r*mie_params['Qext'
bsca = math.sqrt(math.pi/2.0)*(n/np.log(sigma))*r*mie_params['Qsca'
babs = math.sqrt(math.pi/2.0)*(n/np.log(sigma))*r*mie_params['Qabs'
g = math.sqrt(math.pi/2.0)*(n/np.log(sigma))*r*mie_params['Qsca']*mi
bext_array.append(bext) # Appending values to array
bsca_array.append(bsca)
babs_array.append(babs)
g_array.append(g)
return bext_array,bsca_array,babs_array,g_array
bext_array,bsca_array,babs_array,g_array = miecoeff(n,sigma,rad_array,wavele
bext = trap(bext_array,rad_array) # Integrating using trapezo
bsca = trap(bsca_array,rad_array) # rule
babs = trap(babs_array,rad_array)
g = trap(g_array,rad_array)/bsca
print('bext = '+str(bext)+' bsca = '+str(bsca)+' babs = '+str(babs)+' g =
bext = 0.0018502090669298417 bsca = 0.0005138744448271059 babs = 0.001336
334622102735 g = 0.35644716992630454
Visualization of mie parameters
Variation of mie efficiencies with size parameter
In [51]: wavelength = 0.5
ref_index = 1.45+0.05j
rmin = 0.005 # Minimum radius
rmax = 5 # Maximum radius
rad_array = np.arange(rmin,rmax,0.001)
qext_array = []
qsca_array = []
qabs_array = []
size_param_array = []
for r in rad_array:
mie_params = mie.MieQ(ref_index,wavelength,r*2,asDict=True)
size_param = (2*math.pi*r)/wavelength
qext_array.append(mie_params['Qext'])
qabs_array.append(mie_params['Qabs'])
qsca_array.append(mie_params['Qsca'])
size_param_array.append(size_param)
fig = plt.Figure(figsize=(5,5))
plt.plot(size_param_array,qext_array,label = 'Qext')
plt.plot(size_param_array,qabs_array,label = 'Qabs')
plt.plot(size_param_array,qsca_array,label = 'Qsca')
plt.xlabel('Size parameter')
plt.ylabel('Qext/Qabs/Qsca')
plt.legend()
<matplotlib.legend.Legend at 0x7f0d53ff84c0>
Out[51]:
Variation of SSA with refractive index
In [47]: wavelength = 0.5
ref_index = 1.33+0.0001j
rmin = 0.005 # Minimum radius
rmax = 10 # Maximum radius
rad_array = np.arange(rmin,rmax,0.001)
qext_array = []
qsca_array = []
qabs_array = []
size_param_array = []
ssa_array = []
for r in rad_array:
mie_params = mie.MieQ(ref_index,wavelength,r*2,asDict=True)
size_param = (2*math.pi*r)/wavelength
qext_array.append(mie_params['Qext'])
qabs_array.append(mie_params['Qabs'])
qsca_array.append(mie_params['Qsca'])
size_param_array.append(size_param)
ssa_array = np.divide(qsca_array,qext_array)
fig = plt.Figure(figsize=(5,5))
plt.plot(size_param_array,ssa_array,label = str(ref_index))
ref_index = 1.33+0.01j
qext_array = []
qsca_array = []
qabs_array = []
size_param_array = []
ssa_array = []
for r in rad_array:
mie_params = mie.MieQ(ref_index,wavelength,r*2,asDict=True)
size_param = (2*math.pi*r)/wavelength
qext_array.append(mie_params['Qext'])
qabs_array.append(mie_params['Qabs'])
qsca_array.append(mie_params['Qsca'])
size_param_array.append(size_param)
ssa_array = np.divide(qsca_array,qext_array)
plt.plot(size_param_array,ssa_array,label = str(ref_index))
ref_index = 1.33+0.1j
qext_array = []
qsca_array = []
qabs_array = []
size_param_array = []
ssa_array = []
for r in rad_array:
mie_params = mie.MieQ(ref_index,wavelength,r*2,asDict=True)
size_param = (2*math.pi*r)/wavelength
qext_array.append(mie_params['Qext'])
qabs_array.append(mie_params['Qabs'])
qsca_array.append(mie_params['Qsca'])
size_param_array.append(size_param)
ssa_array = np.divide(qsca_array,qext_array)
plt.plot(size_param_array,ssa_array,label = str(ref_index))
plt.xlabel('Size parameter')
plt.ylabel('SSA')
plt.legend()
<matplotlib.legend.Legend at 0x7f0d5416f6a0>
Out[47]:
Useful links
https://scattport.org/index.php/programs-menu/mie-type-codes-menu
http://www.philiplaven.com/mieplot.htm
In [ ]:
