Mie Scattering Algorithms¶
Scott Prahl
Jan 2022
This Jupyter notebook shows the formulas used in miepython
. This code is heavily influenced by Wiscomes MIEV0 code as documented in his paper on Mie scattering and his 1979 NCAR and 1996 NCAR publications.
There are a couple of things that set this code apart from other python Mie codes.
Instead of using the built-in special functions from SciPy, the calculation relies on the logarthmic derivative of the Ricatti-Bessel functions. This technique is significantly more accurate.
The code uses special cases for small spheres. This is faster and more accurate
The code works when the index of refraction
m.real
is zero or whenm.imag
is very large (negative).
The code has been tested up to sizes (\(x=2\pi r/\lambda=10000\)).
If miepython is not installed, uncomment the following cell (i.e., delete the #) and run (shift-enter)
[1]:
#!pip install --user miepython
[2]:
import numpy as np
import matplotlib.pyplot as plt
try:
import miepython.miepython as miepython
except ModuleNotFoundError:
print('miepython not installed. To install, uncomment and run the cell above.')
print('Once installation is successful, rerun this cell again.')
The logarithmic derivative \(D_n\).¶
This routine uses a continued fraction method to compute \(D_n(z)\) proposed by Lentz. Lentz uses the notation :math:`A_n` instead of :math:`D_n`, but I prefer the notation used by Bohren and Huffman. This method eliminates many weaknesses in previous algorithms using forward recursion.
The logarithmic derivative \(D_n\) is defined as
Equation (5) in Lentz’s paper can be used to obtain
Now if
we seek to create
since Lentz showed that
The whole goal is to iterate until the \(\alpha\) and \(\beta\) are identical to the number of digits desired. Once this is achieved, then use equations this equation and the first equation for the logarithmic derivative to calculate \(D_n(z)\).
First terms¶
The value of \(a_j\) is
The first terms for \(\alpha\) and \(\beta\) are then
Later terms¶
To calculate the next \(\alpha\) and \(\beta\), I use
to find the next \(a_j\) and
Calculating \(D_n\)¶
Use formula 7 from Wiscombe’s paper to figure out if upwards or downwards recurrence should be used. Namely if
the upward recurrence would be stable.
The returned array D
is set-up so that \(D_n(z)=\) D[n]
. Therefore the first value for \(D_1(z)\) will not be D[0]
, but rather D[1]
.
\(D_n\) by downwards recurrence.¶
Start downwards recurrence using by accurately calculating D[nstop]
using the Lentz method, then find earlier terms of the logarithmic derivative \(D_n(z)\) using the recurrence relation,
This is a pretty straightforward procedure.
\(D_n\) by upward recurrence.¶
Calculating the logarithmic derivative \(D_n(\rho)\) using the upward recurrence relation,
To calculate the initial value D[1]
we use Wiscombe’s representation that avoids overflow errors when the usual \(D_0(x)=1/tan(z)\) is used.
[3]:
m = 1
x = 1
nstop = 10
dn = np.zeros(nstop, dtype=np.complex128)
print("both techniques work up to 5")
n=5
print(" Lentz",n,miepython._Lentz_Dn(m*x,n).real)
miepython._D_downwards(m*x,nstop, dn)
print("downwards",n, dn[n].real)
miepython._D_upwards(m*x,nstop, dn)
print(" upwards",n, dn[n].real)
print("but upwards fails badly by n=9")
n=9
print(" Lentz",n,miepython._Lentz_Dn(m*x,n).real)
miepython._D_downwards(m*x, nstop, dn)
print("downwards",n,dn[n].real)
miepython._D_upwards(m*x, nstop, dn)
print(" upwards",n,dn[n].real)
both techniques work up to 5
Lentz 5 5.922678838321971
downwards 5 5.922678838321968
upwards 5 5.922678842006903
but upwards fails badly by n=9
Lentz 9 9.952281984945753
downwards 9 9.952281984945756
upwards 9 67.02345742657965
Calculating \(A_n\) and \(B_n\)¶
OK, Here we go. We need to start up the arrays. First, recall (page 128 Bohren and Huffman) that
where \(j_n\) and \(y_n\) are spherical Bessel functions. The first few terms may be worked out as,
and
The main equations for \(a_n\) and \(b_n\) in Bohren and Huffman Equation (4.88).
and
The recurrence relations for \(\psi\) and \(\xi\) depend on the recursion relations for the spherical Bessel functions (page 96 equation 4.11)
where \(z_n\) might be either \(j_n\) or \(y_n\). Thus
If the spheres are perfectly reflecting m.real=0
then Kerker gives equations for \(a_n\) and \(b_n\) that do not depend on \(D_n\) at all
and
Therefore D[n]
will directly correspond to \(D_n\) in Bohren. However, a
and b
will be zero based arrays and so \(a_1\)=a[0]
or \(b_n\)=b[n-1]
[4]:
m=4/3
x=50
print("m=4/3 test, m=",m, " x=",x)
a, b = miepython._mie_An_Bn(m,x)
print("a_1=", a[0])
print("a_1= (0.531105889295-0.499031485631j) #test")
print("b_1=", b[0])
print("b_1= (0.791924475935-0.405931152229j) #test")
print()
m=3/2-1j
x=2
print("upward recurrence test, m=",m, " x=",x)
a, b = miepython._mie_An_Bn(m,x)
print("a_1=", a[0])
print("a_1= (0.546520203397-0.152373857258j) #test")
print("b_1=", b[0])
print("b_1= (0.389714727888+0.227896075256j) #test")
print()
m=11/10-25j
x=2
print("downward recurrence test, m=",m, " x=",x)
a, b = miepython._mie_An_Bn(m,x)
print("a_1=", a[0])
print("a_1= (0.322406907480-0.465063542971j) #test")
print("b_1=", b[0])
print("b_1= (0.575167279092+0.492912495262j) #test")
m=4/3 test, m= 1.3333333333333333 x= 50
a_1= (0.5311058892948326+0.4990314856310949j)
a_1= (0.531105889295-0.499031485631j) #test
b_1= (0.7919244759351927+0.40593115222899945j)
b_1= (0.791924475935-0.405931152229j) #test
upward recurrence test, m= (1.5-1j) x= 2
a_1= (0.5465202033970914+0.1523738572575972j)
a_1= (0.546520203397-0.152373857258j) #test
b_1= (0.3897147278879423-0.22789607525649083j)
b_1= (0.389714727888+0.227896075256j) #test
downward recurrence test, m= (1.1-25j) x= 2
a_1= (0.32240690748075784+0.4650635429715796j)
a_1= (0.322406907480-0.465063542971j) #test
b_1= (0.5751672790921928-0.4929124952616458j)
b_1= (0.575167279092+0.492912495262j) #test
Small Spheres¶
and says this routine should be accurate to six places.
The formula for \({\hat a}_1\) is
where
Note that I have disabled the case when the sphere has no index of refraction. The perfectly conducting sphere equations are
The formula for \({\hat b}_1\) is
The formula for \({\hat a}_2\) is
The scattering and extinction efficiencies are given by
and
with
and the anisotropy (average cosine of the phase function) is
The backscattering efficiency \(Q_\mathrm{back}\) is
where \(S_1(\mu)\) is
[5]:
m=1.5-0.1j
x=0.0665
print("abs(m*x)=",abs(m*x))
qext, qsca, qback, g = miepython._small_mie(m,x)
print("Qext=",qext)
print("Qsca=",qsca)
print("Qabs=",qext-qsca)
print("Qback=",qback)
print("g=",g)
print()
print('The following should be nearly the same as those above:')
print()
x=0.067
print("abs(m*x)=",abs(m*x))
qext, qsca, qback, g = miepython.mie(m,x)
print("Qext=",qext)
print("Qsca=",qsca)
print("Qabs=",qext-qsca)
print("Qback=",qback)
print("g=",g)
abs(m*x)= 0.09997142091617985
Qext= 0.013287673090500258
Qsca= 4.699313232220918e-06
Qabs= 0.013282973777268036
Qback= 7.034282738345809e-06
g= 0.0008751930053081949
The following should be nearly the same as those above:
abs(m*x)= 0.1007230857350985
Qext= 0.013388187040023184
Qsca= 4.842269259059647e-06
Qabs= 0.013383344770764125
Qback= 7.24804468970277e-06
g= 0.0008883995075494251
Small Perfectly Reflecting Spheres¶
The above equations fail when m.real=0
so use these approximations when the sphere is small and refective
[6]:
m = 0 - 0.01j
x=0.099
qext, qsca, qback, g = miepython._small_conducting_mie(m,x)
print("Qext =",qext)
print("Qsca =",qsca)
print("Qabs =",qext-qsca)
print("Qback=",qback)
print("g =",g)
print()
print('The following should be nearly the same as those above:')
print()
m = 0 - 0.01j
x=0.1001
qext, qsca, qback2, g = miepython.mie(m,x)
print("Qext =",qext)
print("Qsca =",qsca)
print("Qabs =",qext-qsca)
print("Qback=",qback2)
print("g =",g)
Qext = 0.0003209674075568899
Qsca = 0.0003209674075568899
Qabs = 0.0
Qback= 0.0008630007227881902
g = -0.3973569106511184
The following should be nearly the same as those above:
Qext = 0.000335472382749474
Qsca = 0.000335472382749474
Qabs = 0.0
Qback= 0.0009019320002161998
g = -0.3973104926046963
Mie scattering calculations¶
From page 120 of Bohren and Huffman the anisotropy is given by
For computation purposes, this must be rewritten as
From page 122 we find an expression for the backscattering efficiency
From page 103 we find an expression for the scattering cross section
The total extinction efficiency is also found on page 103
[7]:
qext, qsca, qback, g = miepython.mie(1.55-0.0j,2*np.pi/0.6328*0.525)
print("Qext=",qext)
print("Qsca=",qsca)
print("Qabs=",qext-qsca)
print("Qback=",qback)
print("g=",g)
Qext= 3.1054255314658765
Qsca= 3.1054255314658765
Qabs= 0.0
Qback= 2.925340649659009
g= 0.6331367580408949
[8]:
x=1000.0
m=1.5-0.1j
qext, qsca, qback, g = miepython.mie(m,x)
print("Qext=",qext)
print("Qsca=",qsca)
print("Qabs=",qext-qsca)
print("Qback=",qback)
print("g=",g)
Qext= 2.0197025208225514
Qsca= 1.1069323889254012
Qabs= 0.9127701318971502
Qback= 0.04153355464458776
g= 0.9508799127402501
[9]:
x=10000.0
m=1.5-1j
qext, qsca, qback, g = miepython.mie(m,x)
print("Qext=",qext)
print("Qsca=",qsca)
print("Qabs=",qext-qsca)
print("Qback=",qback)
print("g=",g)
Qext= 2.004367709696921
Qsca= 1.2365743120721584
Qabs= 0.7677933976247626
Qback= 0.172413800511328
g= 0.8463099581094649
Scattering Matrix¶
The scattering matrix is given by Equation 4.74 in Bohren and Huffman. Namely,
and
If \(\mu=\cos\theta\) then
and
This means that for each angle \(\mu\) we need to know \(\tau_n(\mu)\) and \(\pi_n(\mu)\) for every \(a_n\) and \(b_n\). Equation 4.47 in Bohren and Huffman states
and knowning that \(\pi_0(\mu)=0\) and \(\pi_1(\mu)=1\), all the rest can be found. Similarly
so the plan is to use these recurrence relations to find \(\pi_n(\mu)\) and \(\tau_n(\mu)\) during the summation process.
The only real trick is to account for 0-based arrays when the sums above are 1-based.
[10]:
m=1.55-0.1j
x=5.213
mu = np.array([0.0,0.5,1.0])
S1,S2 = miepython.mie_S1_S2(m,x,mu)
for i in range(len(mu)):
print(mu[i], S2[i].real, S2[i].imag)
0.0 0.0430827037808541 -0.05982417019769807
0.5 -0.08406917971660786 0.13895030592467703
1.0 1.2438024701591521 -0.19843241121775865
Test to match Bohren’s Sample Calculation¶
[11]:
# Test to match Bohren's Sample Calculation
theta = np.arange(0,181,9)
mu=np.cos(theta*np.pi/180)
S1,S2 = miepython.mie_S1_S2(1.55,5.213,mu)
qext, qsca, qback, g = miepython.mie(m,x)
norm = np.sqrt(qext * x**2 * np.pi)
S1 /= norm
S2 /= norm
S11 = (abs(S2)**2 + abs(S1)**2)/2
S12 = (abs(S2)**2 - abs(S1)**2)/2
S33 = (S2 * S1.conjugate()).real
S34 = (S2 * S1.conjugate()).imag
# the minus in POL=-S12/S11 matches that Bohren
# the minus in front of -S34/S11 does not match Bohren's code!
print("ANGLE S11 POL S33 S34")
for i in range(len(mu)):
print("%5d %10.8f % 10.8f % 10.8f % 10.8f" % (theta[i], S11[i]/S11[0], -S12[i]/S11[i], S33[i]/S11[i], -S34[i]/S11[i]))
ANGLE S11 POL S33 S34
0 1.00000000 -0.00000000 1.00000000 -0.00000000
9 0.78538504 -0.00458392 0.99940039 0.03431985
18 0.35688492 -0.04578478 0.98602789 0.16016480
27 0.07660207 -0.36455096 0.84366465 0.39412251
36 0.03553383 -0.53498510 0.68714053 -0.49155756
45 0.07019023 0.00954907 0.95986338 -0.28030538
54 0.05743887 0.04782061 0.98536582 0.16360740
63 0.02196833 -0.44040631 0.64814202 0.62125213
72 0.01259465 -0.83204714 0.20344385 -0.51605054
81 0.01737702 0.03419635 0.79548556 -0.60500689
90 0.01246407 0.23055334 0.93743853 0.26087192
99 0.00679199 -0.71323431 -0.00732217 0.70088744
108 0.00954281 -0.75617653 -0.03954742 -0.65317154
117 0.00863640 -0.28085850 0.53642012 -0.79584669
126 0.00227521 -0.23864148 0.96777914 0.08033545
135 0.00544047 -0.85116040 0.18710096 -0.49042758
144 0.01602875 -0.70649116 0.49501921 -0.50579267
153 0.01889077 -0.89109951 0.45322894 -0.02291691
162 0.01952522 -0.78348591 -0.39140822 0.48264836
171 0.03016127 -0.19626673 -0.96204724 0.18959028
180 0.03831054 -0.00000000 -1.00000000 -0.00000000
[12]:
num=100
m=1.1
x=np.linspace(0.01,0.21,num)
qext, qsca, qback, g = miepython.mie(m,x)
plt.plot(x,qback)
plt.plot((abs(0.1/m),abs(0.1/m)),(0,qback[num-1]))
plt.xlabel("Size Parameter (-)")
plt.ylabel("Backscattering Efficiency")
plt.show()
[ ]: