In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import miepython as mp
Mie scattering describes the special case of the interaction of light passing through a non-absorbing medium with a single embedded spherical object. The sphere itself can be non-absorbing, moderately absorbing, or perfectly absorbing.
Specifically, the scattering function $p(\theta_i,\phi_i,\theta_o,\phi_o)$ describes the amount of light scattered by a particle for light incident at an angle $(\theta_i,\phi_i)$ and exiting the particle (in the far field) at an angle $(\theta_o,\phi_o)$. For simplicity, the scattering function is often assumed to be rotationally symmetric (it is, obviously, for spherical scatterers) and that the angle that the light is scattered into only depends the $\theta=\theta_o-\theta_i$. In this case, the scattering function can be written as $p(\theta)$. Finally, the angle is often replaced by $\mu=\cos\theta$ and therefore the phase function becomes just $p(\mu)$.
The figure below shows the basic idea. An incoming monochromatic plane wave hits a sphere and produces in the far field two separate monochromatic waves — a slightly attenuated unscattered planar wave and an outgoing spherical wave.
Obviously. the scattered light will be cylindrically symmetric about the ray passing through the center of the sphere.
In [2]:
t = np.linspace(0,2*np.pi,100)
xx = np.cos(t)
yy = np.sin(t)
fig,ax=plt.subplots(figsize=(10,8))
plt.axes().set_aspect('equal')
plt.plot(xx,yy)
plt.plot([-5,7],[0,0],'--k')
plt.annotate('incoming irradiance', xy=(-4.5,-2.3),ha='left',color='blue',fontsize=14)
for i in range(6):
y0 = i -2.5
plt.annotate('',xy=(-1.5,y0),xytext=(-5,y0),arrowprops=dict(arrowstyle="->",color='blue'))
plt.annotate('unscattered irradiance', xy=(3,-2.3),ha='left',color='blue',fontsize=14)
for i in range(6):
y0 = i -2.5
plt.annotate('',xy=(7,y0),xytext=(3,y0),arrowprops=dict(arrowstyle="->",color='blue',ls=':'))
plt.annotate('scattered\nspherical\nwave', xy=(0,1.5),ha='left',color='red',fontsize=16)
plt.annotate('',xy=(2.5,2.5),xytext=(0,0),arrowprops=dict(arrowstyle="->",color='red'))
plt.annotate(r'$\theta$',xy=(2,0.7),color='red',fontsize=14)
plt.annotate('',xy=(2,2),xytext=(2.7,0),arrowprops=dict(connectionstyle="arc3,rad=0.2", arrowstyle="<->",color='red'))
plt.xlim(-5,7)
plt.ylim(-3,3)
plt.axis('off')
plt.show()
In [3]:
fig,ax=plt.subplots(figsize=(10,8))
plt.axes().set_aspect('equal')
plt.scatter([0],[0],s=30)
m = 1.5
x = np.pi/3
theta = np.linspace(-180,180,180)
mu = np.cos(theta/180*np.pi)
scat = 15 * mp.i_unpolarized(m,x,mu)
plt.plot(scat*np.cos(theta/180*np.pi),scat*np.sin(theta/180*np.pi))
for i in range(12):
ii = i*15
xx = scat[ii]*np.cos(theta[ii]/180*np.pi)
yy = scat[ii]*np.sin(theta[ii]/180*np.pi)
# print(xx,yy)
plt.annotate('',xy=(xx,yy),xytext=(0,0),arrowprops=dict(arrowstyle="->",color='red'))
plt.annotate('incident irradiance', xy=(-4.5,-2.3),ha='left',color='blue',fontsize=14)
for i in range(6):
y0 = i -2.5
plt.annotate('',xy=(-1.5,y0),xytext=(-5,y0),arrowprops=dict(arrowstyle="->",color='blue'))
plt.annotate('unscattered irradiance', xy=(3,-2.3),ha='left',color='blue',fontsize=14)
for i in range(6):
y0 = i -2.5
plt.annotate('',xy=(7,y0),xytext=(3,y0),arrowprops=dict(arrowstyle="->",color='blue',ls=':'))
plt.annotate('scattered\nspherical wave', xy=(0,1.5),ha='left',color='red',fontsize=16)
plt.xlim(-5,7)
plt.ylim(-3,3)
#plt.axis('off')
plt.show()
So the scattering function or phase function has at least three reasonable normalizations that involve integrating over all $4\pi$ steradians. Below $d\Omega=\sin\theta d\theta\,d\phi$ is a differential solid angle
$$ \begin{align} \int_{4\pi} p(\theta,\phi) \,d\Omega &= 1\\[2mm] \int_{4\pi} p(\theta,\phi) \,d\Omega &= 4\pi \\[2mm] \int_{4\pi} p(\theta,\phi) \,d\Omega &= a \qquad\qquad \mbox{Used by miepython}\\[2mm] \end{align} $$where $a$ is the single scattering albedo,
$$ a = \frac{\sigma_s}{\sigma_s+\sigma_a} $$and $\sigma_s$ is the scattering cross section, and $\sigma_a$ is the absorption cross section.
The choice of normalization was made because it accounts for light lost through absorption by the sphere.
If the incident light has units of watts, then the values from the scattering function $p(\theta,\phi)$ have units of radiant intensity or W/sr.
For example, a circular detector with radius $r_d$ at a distance $R$ will subtend an angle
$$ \Omega = \frac{\pi r_d^2}{R^2} $$(assuming $r_d\ll R$). Now if $P_0$ of light is scattered by a sphere then the scattered power on the detector will be
$$ P_d = P_0 \cdot \Omega \cdot p(\theta,\phi) $$
In [4]:
m = 1.5
x = np.pi/3
theta = np.linspace(-180,180,180)
mu = np.cos(theta/180*np.pi)
scat = mp.i_unpolarized(m,x,mu)
fig,ax = plt.subplots(1,2,figsize=(12,5))
ax=plt.subplot(121, projection='polar')
ax.plot(theta/180*np.pi,scat)
ax.set_rticks([0.05, 0.1,0.15])
ax.set_title("m=1.5, Sphere Diameter = $\lambda$/3")
plt.subplot(122)
plt.plot(theta,scat)
plt.xlabel('Exit Angle [degrees]')
plt.ylabel('Unpolarized Scattered light [1/sr]')
plt.title('m=1.5, Sphere Diameter = $\lambda$/3')
plt.ylim(0.00,0.2)
plt.show()
Classic Rayleigh scattering treats small particles with natural (unpolarized) light.
The solid black line denotes the total scattered intensity. The red dashed line is light scattered that is polarized perpendicular to the plane of the graph and the blue dotted line is for light parallel to the plane of the graph. (Compare with van de Hult, Figure 10)
In [5]:
m = 1.3
x = 0.01
theta = np.linspace(-180,180,180)
mu = np.cos(theta/180*np.pi)
ipar = mp.i_par(m,x,mu)/2
iper = mp.i_per(m,x,mu)/2
iun = mp.i_unpolarized(m,x,mu)
fig,ax = plt.subplots(1,2,figsize=(12,5))
ax=plt.subplot(121, projection='polar')
ax.plot(theta/180*np.pi,iper,'r--')
ax.plot(theta/180*np.pi,ipar,'b:')
ax.plot(theta/180*np.pi,iun,'k')
ax.set_rticks([0.05, 0.1,0.15])
plt.title('m=%.2f, Sphere Parameter = %.2f' %(m,x))
plt.subplot(122)
plt.plot(theta,iper,'r--')
plt.plot(theta,ipar,'b:')
plt.plot(theta,iun,'k')
plt.xlabel('Exit Angle [degrees]')
plt.ylabel('Normalized Scattered light [1/sr]')
plt.title('m=%.2f, Sphere Parameter = %.2f' %(m,x))
plt.ylim(0.00,0.125)
plt.text(130,0.02,r"$0.5I_{per}$",color="blue", fontsize=16)
plt.text(120,0.062,r"$0.5I_{par}$",color="red", fontsize=16)
plt.text(30,0.11,r"$I_{unpolarized}$",color="black", fontsize=16)
plt.show()
Sometimes one would like the scattering function normalized so that the integral over all $4\pi$ steradians to be the scattering cross section
$$ \sigma_{sca} = \frac{\pi d^2}{4} Q_{sca} $$The differential scattering cross section \frac{d\sigma_{sca}}{d\Omega}
$$ \sigma_{sca} = \int_{4\pi} \frac{d\sigma_{sca}}{d\Omega}\,d\Omega $$Since the unpolarized scattering normalized so its integral is the single scattering albedo, this means that
$$ \frac{Q_{sca}}{Q_{ext}} = \int_{4\pi} p(\mu) \sin\theta\,d\theta d\phi $$and therefore the differential scattering cross section can be obtained miepython
using
Note that this is $Q_{ext}$ and not $Q_{sca}$ because of the choice of normalization!
For example, here is a replica of figure 4
In [6]:
m = 1.4-0j
lambda0 = 532e-9 # m
theta = np.linspace(0,180,1000)
mu = np.cos(theta* np.pi/180)
d = 1700e-9 # m
x = 2 * np.pi/lambda0 * d/2
geometric_cross_section = np.pi * d**2/4 * 1e4 # cm**2
qext, qsca, qback, g = mp.mie(m,x)
sigma_sca = geometric_cross_section * qext * mp.i_unpolarized(m,x,mu)
plt.semilogy(theta, sigma_sca*1e-3, color='blue')
plt.text(15, sigma_sca[0]*3e-4, "%.0fnm\n(x10$^{-3}$)" % (d*1e9), color='blue')
d = 170e-9 # m
x = 2 * np.pi/lambda0 * d/2
geometric_cross_section = np.pi * d**2/4 * 1e4 # cm**2
qext, qsca, qback, g = mp.mie(m,x)
sigma_sca = geometric_cross_section * qext * mp.i_unpolarized(m,x,mu)
plt.semilogy(theta, sigma_sca, color='red')
plt.text(110, sigma_sca[-1]/2, "%.0fnm" % (d*1e9), color='red')
d = 17e-9 # m
x = 2 * np.pi/lambda0 * d/2
geometric_cross_section = np.pi * d**2/4 * 1e4 # cm**2
qext, qsca, qback, g = mp.mie(m,x)
sigma_sca = geometric_cross_section * qext * mp.i_unpolarized(m,x,mu)
plt.semilogy(theta, sigma_sca*1e6, color='green')
plt.text(130, sigma_sca[-1]*1e6, "(x10$^6$)\n%.0fnm" % (d*1e9), color='green')
plt.title("Refractive index m=1.4, $\lambda$=532nm")
plt.xlabel("Scattering Angle (degrees)")
plt.ylabel("Diff. Scattering Cross Section (cm$^2$/sr)")
plt.grid(True)
plt.show()
Specifically, to ensure proper normalization, the integral of the scattering function over all solid angles must be unity
$$ a = \int_0^{2\pi}\int_0^\pi \, p(\theta,\phi)\,\sin\theta\,d\theta\,d\phi $$or with a change of variables $\mu=\cos\theta$ and using the symmetry to the integral in $\phi$
$$ a = 2\pi \int_{-1}^1 \, p(\mu)\,d\mu $$This integral can be done numerically by simply summing all the rectangles
$$ a = 2\pi \sum_{i=0}^N p(\mu_i)\,\Delta\mu_i $$and if all the rectanges have the same width
$$ a = 2\pi\Delta\mu \sum_{i=0}^N p(\mu_i) $$
In [7]:
m = 1.5
x = 1
mu = np.linspace(-1,1,501)
intensity = mp.i_unpolarized(m,x,mu)
qext, qsca, qback, g = mp.mie(m,x)
a = qsca/qext
#integrate over all angles
dmu = mu[1] - mu[0]
total = 2 * np.pi * dmu * np.sum(intensity)
plt.plot(mu,intensity)
plt.xlabel(r'$\cos(\theta)$')
plt.ylabel('Unpolarized Scattering Intensity [1/sr]')
plt.title('m=%.3f%+.3fj, x=%.2f, a=%.3f, total=%.3f'%(m.real,m.imag,x,a, total))
plt.show()
In [8]:
m = 1.5 - 1.5j
x = 1
mu = np.linspace(-1,1,501)
intensity = mp.i_unpolarized(m,x,mu)
qext, qsca, qback, g = mp.mie(m,x)
a = qsca/qext
#integrate over all angles
dmu = mu[1] - mu[0]
total = 2 * np.pi * dmu * np.sum(intensity)
plt.plot(mu,intensity)
plt.xlabel(r'$\cos(\theta)$')
plt.ylabel('Unpolarized Scattering Intensity [1/sr]')
plt.title('m=%.3f%+.3fj, x=%.2f, a=%.3f, total=%.3f'%(m.real,m.imag,x,a, total))
plt.show()
In [9]:
m = 1.5-1.5j
x = 1
theta = np.linspace(0,180,361)*np.pi/180
mu = np.cos(theta)
intensity = mp.i_unpolarized(m,x,mu)
qext, qsca, qback, g = mp.mie(m,x)
a = qsca/qext
#integrate over all angles
dtheta = theta[1]-theta[0]
total = 2 * np.pi * dtheta * np.sum(intensity* np.sin(theta))
plt.plot(mu,intensity)
plt.xlabel(r'$\cos(\theta)$')
plt.ylabel('Unpolarized Scattering Intensity [1/sr]')
plt.title('m=%.3f%+.3fj, x=%.2f, a=%.3f, total=%.3f'%(m.real,m.imag,x,a, total))
plt.show()
Wiscombe normalizes as
$$ \int_{4\pi} p(\theta,\phi) \,d\Omega = \pi x^2 Q_{sca} $$where $p(\theta)$ is the scattered light.
Once corrected for differences in phase function normalization, Wiscombe's test cases match those from miepython
exactly.
In [10]:
"""
MIEV0 Test Case 14: Refractive index: real 1.500 imag -1.000E+00, Mie size parameter = 1.000
Angle Cosine S-sub-1 S-sub-2 Intensity Deg of Polzn
0.00 1.000000 5.84080E-01 1.90515E-01 5.84080E-01 1.90515E-01 3.77446E-01 0.0000
30.00 0.866025 5.65702E-01 1.87200E-01 5.00161E-01 1.45611E-01 3.13213E-01 -0.1336
60.00 0.500000 5.17525E-01 1.78443E-01 2.87964E-01 4.10540E-02 1.92141E-01 -0.5597
90.00 0.000000 4.56340E-01 1.67167E-01 3.62285E-02 -6.18265E-02 1.20663E-01 -0.9574
"""
x=1.0
m=1.5-1.0j
mu=np.cos(np.linspace(0,90,4) * np.pi/180)
qext, qsca, qback, g = mp.mie(m,x)
albedo = qsca/qext
unpolar = mp.i_unpolarized(m,x,mu) # normalized to a
unpolar /= albedo # normalized to 1
unpolar_miev = np.array([3.77446E-01,3.13213E-01,1.92141E-01,1.20663E-01])
unpolar_miev /= np.pi * qsca * x**2 # normalized to 1
ratio = unpolar_miev/unpolar
print("MIEV0 Test Case 14: m=1.500-1.000j, Mie size parameter = 1.000")
print()
print(" %9.1f°%9.1f°%9.1f°%9.1f°"%(0,30,60,90))
print("MIEV0 %9.5f %9.5f %9.5f %9.5f"%(unpolar_miev[0],unpolar_miev[1],unpolar_miev[2],unpolar_miev[3]))
print("miepython %9.5f %9.5f %9.5f %9.5f"%(unpolar[0],unpolar[1],unpolar[2],unpolar[3]))
print("ratio %9.5f %9.5f %9.5f %9.5f"%(ratio[0],ratio[1],ratio[2],ratio[3]))
In [11]:
"""
MIEV0 Test Case 10: Refractive index: real 1.330 imag -1.000E-05, Mie size parameter = 100.000
Angle Cosine S-sub-1 S-sub-2 Intensity Deg of Polzn
0.00 1.000000 5.25330E+03 -1.24319E+02 5.25330E+03 -1.24319E+02 2.76126E+07 0.0000
30.00 0.866025 -5.53457E+01 -2.97188E+01 -8.46720E+01 -1.99947E+01 5.75775E+03 0.3146
60.00 0.500000 1.71049E+01 -1.52010E+01 3.31076E+01 -2.70979E+00 8.13553E+02 0.3563
90.00 0.000000 -3.65576E+00 8.76986E+00 -6.55051E+00 -4.67537E+00 7.75217E+01 -0.1645
"""
x=100.0
m=1.33-1e-5j
mu=np.cos(np.linspace(0,90,4) * np.pi/180)
qext, qsca, qback, g = mp.mie(m,x)
albedo = qsca/qext
unpolar = mp.i_unpolarized(m,x,mu) # normalized to a
unpolar /= albedo # normalized to 1
unpolar_miev = np.array([2.76126E+07,5.75775E+03,8.13553E+02,7.75217E+01])
unpolar_miev /= np.pi * qsca * x**2 # normalized to 1
ratio = unpolar_miev/unpolar
print("MIEV0 Test Case 10: m=1.330-0.00001j, Mie size parameter = 100.000")
print()
print(" %9.1f°%9.1f°%9.1f°%9.1f°"%(0,30,60,90))
print("MIEV0 %9.5f %9.5f %9.5f %9.5f"%(unpolar_miev[0],unpolar_miev[1],unpolar_miev[2],unpolar_miev[3]))
print("miepython %9.5f %9.5f %9.5f %9.5f"%(unpolar[0],unpolar[1],unpolar[2],unpolar[3]))
print("ratio %9.5f %9.5f %9.5f %9.5f"%(ratio[0],ratio[1],ratio[2],ratio[3]))
In [12]:
"""
MIEV0 Test Case 7: Refractive index: real 0.750 imag 0.000E+00, Mie size parameter = 10.000
Angle Cosine S-sub-1 S-sub-2 Intensity Deg of Polzn
0.00 1.000000 5.58066E+01 -9.75810E+00 5.58066E+01 -9.75810E+00 3.20960E+03 0.0000
30.00 0.866025 -7.67288E+00 1.08732E+01 -1.09292E+01 9.62967E+00 1.94639E+02 0.0901
60.00 0.500000 3.58789E+00 -1.75618E+00 3.42741E+00 8.08269E-02 1.38554E+01 -0.1517
90.00 0.000000 -1.78590E+00 -5.23283E-02 -5.14875E-01 -7.02729E-01 1.97556E+00 -0.6158
"""
x=10.0
m=0.75
mu=np.cos(np.linspace(0,90,4) * np.pi/180)
qext, qsca, qback, g = mp.mie(m,x)
albedo = qsca/qext
unpolar = mp.i_unpolarized(m,x,mu) # normalized to a
unpolar /= albedo # normalized to 1
unpolar_miev = np.array([3.20960E+03,1.94639E+02,1.38554E+01,1.97556E+00])
unpolar_miev /= np.pi * qsca * x**2 # normalized to 1
ratio = unpolar_miev/unpolar
print("MIEV0 Test Case 7: m=0.75, Mie size parameter = 10.000")
print()
print(" %9.1f°%9.1f°%9.1f°%9.1f°"%(0,30,60,90))
print("MIEV0 %9.5f %9.5f %9.5f %9.5f"%(unpolar_miev[0],unpolar_miev[1],unpolar_miev[2],unpolar_miev[3]))
print("miepython %9.5f %9.5f %9.5f %9.5f"%(unpolar[0],unpolar[1],unpolar[2],unpolar[3]))
print("ratio %9.5f %9.5f %9.5f %9.5f"%(ratio[0],ratio[1],ratio[2],ratio[3]))
In [13]:
"""
BHMie Test Case 14, Refractive index = 1.5000-1.0000j, Size parameter = 1.0000
Angle Cosine S1 S2
0.00 1.0000 -8.38663e-01 -8.64763e-01 -8.38663e-01 -8.64763e-01
0.52 0.8660 -8.19225e-01 -8.61719e-01 -7.21779e-01 -7.27856e-01
1.05 0.5000 -7.68157e-01 -8.53697e-01 -4.19454e-01 -3.72965e-01
1.57 0.0000 -7.03034e-01 -8.43425e-01 -4.44461e-02 6.94424e-02
"""
x=1.0
m=1.5-1j
mu=np.cos(np.linspace(0,90,4) * np.pi/180)
qext, qsca, qback, g = mp.mie(m,x)
albedo = qsca/qext
unpolar = mp.i_unpolarized(m,x,mu) # normalized to a
unpolar /= albedo # normalized to 1
s1_bh = np.empty(4,dtype=np.complex)
s1_bh[0] = -8.38663e-01 - 8.64763e-01*1j
s1_bh[1] = -8.19225e-01 - 8.61719e-01*1j
s1_bh[2] = -7.68157e-01 - 8.53697e-01*1j
s1_bh[3] = -7.03034e-01 - 8.43425e-01*1j
s2_bh = np.empty(4,dtype=np.complex)
s2_bh[0] = -8.38663e-01 - 8.64763e-01*1j
s2_bh[1] = -7.21779e-01 - 7.27856e-01*1j
s2_bh[2] = -4.19454e-01 - 3.72965e-01*1j
s2_bh[3] = -4.44461e-02 + 6.94424e-02*1j
# BHMie seems to normalize their intensities to 4 * pi * x**2 * Qsca
unpolar_bh = (abs(s1_bh)**2+abs(s2_bh)**2)/2
unpolar_bh /= np.pi * qsca * 4 * x**2 # normalized to 1
ratio = unpolar_bh/unpolar
print("BHMie Test Case 14: m=1.5000-1.0000j, Size parameter = 1.0000")
print()
print(" %9.1f°%9.1f°%9.1f°%9.1f°"%(0,30,60,90))
print("BHMIE %9.5f %9.5f %9.5f %9.5f"%(unpolar_bh[0],unpolar_bh[1],unpolar_bh[2],unpolar_bh[3]))
print("miepython %9.5f %9.5f %9.5f %9.5f"%(unpolar[0],unpolar[1],unpolar[2],unpolar[3]))
print("ratio %9.5f %9.5f %9.5f %9.5f"%(ratio[0],ratio[1],ratio[2],ratio[3]))
print()
print("Note that this test is identical to MIEV0 Test Case 14 above.")
print()
print("Wiscombe's code is much more robust than Bohren's so I attribute errors all to Bohren")
Tiny water droplet (0.26 microns) in clouds has pretty strong forward scattering! A graph of this is figure 4.9 in Bohren and Huffman's Absorption and Scattering of Light by Small Particles.
A bizarre scaling factor of $16\pi$ is needed to make the miepython
results match those in the figure 4.9.
In [14]:
x=3
m=1.33-1e-8j
theta = np.linspace(0,180,181)
mu = np.cos(theta*np.pi/180)
scaling_factor = 16*np.pi
iper = scaling_factor*mp.i_per(m,x,mu)
ipar = scaling_factor*mp.i_par(m,x,mu)
P = (iper-ipar)/(iper+ipar)
plt.subplots(2,1,figsize=(8,8))
plt.subplot(2,1,1)
plt.semilogy(theta,ipar,label='$i_{par}$')
plt.semilogy(theta,iper,label='$i_{per}$')
plt.xlim(0,180)
plt.xticks(range(0,181,30))
plt.ylabel('i$_{par}$ and i$_{per}$')
plt.legend()
plt.title('Figure 4.9 from Bohren & Huffman')
plt.subplot(2,1,2)
plt.plot(theta,P)
plt.ylim(-1,1)
plt.xticks(range(0,181,30))
plt.xlim(0,180)
plt.ylabel('Polarization')
plt.plot([0,180],[0,0],':k')
plt.xlabel('Angle (Degrees)')
plt.show()
In [15]:
x=5
m=10000
theta = np.linspace(0,180,361)
mu = np.cos(theta*np.pi/180)
fig, ax = plt.subplots(figsize=(8,8))
x=10
s1,s2 = mp.mie_S1_S2(m,x,mu)
sone = 2.5*abs(s1)
stwo = 2.5*abs(s2)
plt.plot(theta,sone,'b')
plt.plot(theta,stwo,'--r')
plt.annotate('x=%.1f '%x,xy=(theta[-1],sone[-1]),ha='right',va='bottom')
x=5
s1,s2 = mp.mie_S1_S2(m,x,mu)
sone = 2.5*abs(s1) + 1
stwo = 2.5*abs(s2) + 1
plt.plot(theta,sone,'b')
plt.plot(theta,stwo,'--r')
plt.annotate('x=%.1f '%x,xy=(theta[-1],sone[-1]),ha='right',va='bottom')
x=3
s1,s2 = mp.mie_S1_S2(m,x,mu)
sone = 2.5*abs(s1) + 2
stwo = 2.5*abs(s2) + 2
plt.plot(theta,sone,'b')
plt.plot(theta,stwo,'--r')
plt.annotate('x=%.1f '%x,xy=(theta[-1],sone[-1]),ha='right',va='bottom')
x=1
s1,s2 = mp.mie_S1_S2(m,x,mu)
sone = 2.5*abs(s1) + 3
stwo = 2.5*abs(s2) + 3
plt.plot(theta,sone,'b')
plt.plot(theta,stwo,'--r')
plt.annotate('x=%.1f '%x,xy=(theta[-1],sone[-1]),ha='right',va='bottom')
x=0.5
s1,s2 = mp.mie_S1_S2(m,x,mu)
sone = 2.5*abs(s1) + 4
stwo = 2.5*abs(s2) + 4
plt.plot(theta,sone,'b')
plt.plot(theta,stwo,'--r')
plt.annotate('x=%.1f '%x,xy=(theta[-1],sone[-1]),ha='right',va='bottom')
plt.xlim(0,180)
plt.ylim(0,5.5)
plt.xticks(range(0,181,30))
plt.yticks(np.arange(0,5.51,0.5))
plt.title('Figure 29 from van de Hulst, Non-Absorbing Spheres')
plt.xlabel('Angle (Degrees)')
ax.set_yticklabels(['0','1/2','0','1/2','0','1/2','0','1/2','0','1/2','5',' '])
plt.grid(True)
plt.show()
In [16]:
## Kerker, Angular Gain
x=1
m=10000
theta = np.linspace(0,180,361)
mu = np.cos(theta*np.pi/180)
fig, ax = plt.subplots(figsize=(8,8))
s1,s2 = mp.mie_S1_S2(m,x,mu)
G1 = 4*abs(s1)**2/x**2
G2 = 4*abs(s2)**2/x**2
plt.plot(theta,G1,'b')
plt.plot(theta,G2,'--r')
plt.annotate('$G_1$',xy=(50,0.36),color='blue',fontsize=14)
plt.annotate('$G_2$',xy=(135,0.46),color='red',fontsize=14)
plt.xlim(0,180)
plt.xticks(range(0,181,30))
plt.title('Figure 4.51 from Kerker, Non-Absorbing Spheres, x=1')
plt.xlabel('Angle (Degrees)')
plt.ylabel('Angular Gain')
plt.show()
In [ ]: