Mie Random Deviates¶
Scott Prahl
Mar 2021, Version 4
The problem is to generate random scattering angles which match a given Mie scattering profile.
This is difficult when the size parameter gets large because nearly all the light is scattered directly forward.
This notebook is an attempt to solve the scattering problem.
If RigolWFM 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
except ModuleNotFoundError:
print('miepython not installed. To install, uncomment and run the cell above.')
print('Once installation is successful, rerun this cell again.')
Random Deviates from a PDF¶
One method of generating a random number \(\xi\) with a specified distribution \(p(\xi)\) is to create a random event for the variable \(\xi\) such that the random event falls with frequency \(p(x)dx\) in the interval \((\xi,\xi+d\xi)\). This method requires the normalization of the probability density function (PDF) over the interval \((a,b)\)
This is done by choosing a random number \(R\) uniformly distributed in the interval \([0,1]\) and requiring
Note that \(R(\xi)\) represents the cumulative distribution function for \(p(\xi')\).
Azimuthal Angles¶
A normalized phase function describes the probability density function for the azimuthal and longitudinal angles for a photon when it is scattered. If the phase function has no azimuthal dependence, then the azimuthal angle \(\phi\) is uniformly distributed between 0 and \(2\pi\), and may be generated by multiplying a pseudo-random number \(R\) uniformly distributed over the interval [0,1] by \(2\pi\)
Uniformly distributed longitudinal angles¶
The probability density function for the longitudinal angle \(\theta\) between the current photon direction and the scattered photon direction is found by integrating the phase function over all azimuthal angles \(p(\cos\theta)\). For example, the probability density function for an isotropic distribution is
Substituting Equation (A1.9) into Equation (A1.2) yields the following generating function for cosine of the longitudinal angle \(\theta\)
Random Deviates for Henyey Greenstein¶
The probability density function corresponding to the Henyey-Greenstein phase function is
This equation should not be used for isotropic scattering — because of division by zero — use the equation above.
Cumulative Distribution Function for Mie scattering¶
Unfortunately, we cannot do the integral and solve for the angle analytically for the Mie scattering case. We will need to do it numerically.
We ignore polarization effects and are just considering total scattering of the electric field. Moreover, we assume that the scattering function \(S(\theta,\phi)\) is independent of azimuthal angle \(\phi\). In that case the cumulative distribution function (CDF) is
if \(\mu=\cos\theta\). Then
and of course \(\mathrm{CDF}(-1)=0\) and \(\mathrm{CDF}(1)=1\).
The unpolarized scattering function¶
Consider relatively isotropic mie scattering from a smallish sphere. The idea is that we want to generate random angles with this distribution.
[3]:
lambdaa = 0.550 # microns
m = 1.33 # water
r = 0.15 # microns
x = 2*np.pi*r/lambdaa
num=50
mu = np.linspace(-1,1,num)
s1, s2 = miepython.mie_S1_S2(m,x,mu)
scat = (abs(s1)**2 + abs(s2)**2)/2
plt.scatter(mu, scat,s=1)
plt.xlabel('Cosine of Exit Angle')
plt.ylabel('Unpolarized Scattering Function')
plt.title(r'Water Droplet ($\lambda$=550nm, r=%.2f$\mu$m)'%r)
plt.ylim([-0.1,0.4])
plt.xlim([-1.1,1.1])
plt.show()
The CDF¶
miepython has a function to generate the CDF directly
[4]:
help(miepython.mie_cdf)
Help on function mie_cdf in module miepython.miepython:
mie_cdf(m, x, num)
Create a CDF for unpolarized scattering uniformly spaced in cos(theta).
The CDF covers scattered (exit) angles ranging from 180 to 0 degrees.
(The cosines are uniformly distributed over -1 to 1.) Because the angles
are uniformly distributed in cos(theta), the scattering function is not
sampled uniformly and therefore huge array sizes are needed to adequately
sample highly anisotropic phase functions.
Since this is a cumulative distribution function, the maximum value
should be 1.
Args:
m: the complex index of refraction of the sphere
x: the size parameter of the sphere
num: length of desired CDF array
Returns:
mu: array of cosines of angles
cdf: array of cumulative distribution function values
[5]:
lambdaa = 0.550 # microns
m = 1.33 # water
r = 0.15 # microns
x = 2*np.pi*r/lambdaa
num=50
mu, cdf = miepython.mie_cdf(m,x,num)
plt.scatter(mu,cdf,s=1)
plt.xlabel('Cosine of Exit Angle')
plt.ylabel('CDF of Unpolarized Scattering Function')
plt.title(r'Water Droplet ($\lambda$=550nm, r=%.2f$\mu$m)'%r)
plt.ylim([-0.1,1.1])
plt.xlim([-1.1,1.1])
plt.show()
Inverting¶
To solve, we just reverse the \(x\) and \(y\) axes.
[6]:
lambdaa = 0.550 # microns
m = 1.33 # water
r = 0.15 # microns
x = 2*np.pi*r/lambdaa
num=20
mu, cdf = miepython.mie_cdf(m,x,num)
plt.scatter(cdf, mu, s=1)
plt.ylabel('Cosine of Exit Angle')
plt.xlabel('CDF of Unpolarized Scattering Function')
plt.title(r'Water Droplet ($\lambda$=550nm, r=%.2f$\mu$m)'%r)
plt.xlim([-0.1,1.1])
plt.ylim([-1.1,1.1])
plt.show()
Better Inversion¶
In the graph above, the horizontal spacing is not uniform. For speed in the Monte Carlo program we would like direct look ups into the array.
Fortunately there is another function that returns a CDF with uniform spacing in CDF
[7]:
help(miepython.mie_mu_with_uniform_cdf)
Help on function mie_mu_with_uniform_cdf in module miepython.miepython:
mie_mu_with_uniform_cdf(m, x, num)
Create a CDF for unpolarized scattering for uniform CDF.
The CDF covers scattered (exit) angles ranging from 180 to 0 degrees.
(The cosines are uniformly distributed over -1 to 1.) These angles mu
correspond to uniform spacing of the cumulative distribution function
for unpolarized Mie scattering where cdf[i] = i/(num-1).
This is a brute force implementation that solves the problem by
calculating the CDF at many points and then scanning to find the
specific angles that correspond to uniform interval of the CDF.
Since this is a cumulative distribution function, the maximum value
should be 1.
Args:
m: the complex index of refraction of the sphere
x: the size parameter of the sphere
num: length of desired CDF array
Returns:
mu: array of cosines of angles (irregularly spaced)
cdf: array of cumulative distribution function values
[8]:
lambdaa = 0.550 # microns
m = 1.33 # water
r = 0.15 # microns
x = 2*np.pi*r/lambdaa
num=20
mid = num // 2
mu, cdf = miepython.mie_mu_with_uniform_cdf(m,x,num)
plt.scatter(cdf, mu, s=1)
plt.plot([cdf[mid],cdf[mid]],[-1.1,mu[mid]],color='blue')
plt.plot([-0.1,cdf[mid]],[mu[mid],mu[mid]],color='blue')
plt.scatter([cdf[mid],-0.1],[-1.1,mu[mid]],s=20)
plt.ylabel('Cosine of Exit Angle')
plt.xlabel('CDF of Unpolarized Scattering Function')
plt.title(r'Water Droplet ($\lambda$=550nm, r=%.2f$\mu$m)'%r)
plt.xlim([-0.1,1.1])
plt.ylim([-1.1,1.1])
plt.show()
#print(cdf[mid],mu[mid])
And now the a random deviate along (say 0.526) will be mapped to the proper exit angle (0.715)
Generating random Mie scattering angles¶
So now once we have calculated the magic mu array that has a CDF that is uniformly spaced, we can just do a quick look up to get the next random deviate.
[9]:
help(miepython.generate_mie_costheta)
Help on function generate_mie_costheta in module miepython.miepython:
generate_mie_costheta(mu_cdf)
Generate a new scattering angle using a cdf.
A uniformly spaced cumulative distribution function (CDF) is needed.
New random angles are generated by selecting a random interval
mu[i] to mu[i+1] and choosing an angle uniformly distributed over
the interval.
Args:
mu_cdf: a cumulative distribution function
Returns
The cosine of the scattering angle
[10]:
lambdaa = 0.550 # microns
m = 1.33 # water
r = 0.15 # microns
x = 2*np.pi*r/lambdaa
num_angles=20
mu, cdf = miepython.mie_mu_with_uniform_cdf(m,x,num_angles)
# calculate the phase function at each angle
s1,s2 = miepython.mie_S1_S2(m,x,mu)
s = (abs(s1)**2+abs(s2)**2)/2
# generate a bunch of random angles
num_deviates = 100000
angles = np.empty(num_deviates)
for i in range(num_deviates) :
angles[i] = miepython.generate_mie_costheta(mu)
num_bins = 20
plt.hist(angles, bins=num_bins)
plt.plot(mu,s*num_deviates/num_bins*4*np.pi)
#plt.yscale('log')
plt.title("r=%.2f$\mu$m, m=%.3f and %d bins"%(r,m.real,num_bins))
plt.xlabel('Cosine of Exit Angle')
plt.ylabel('Frequency')
plt.xlim([-1.1,1.1])
plt.show()
[ ]: