<p>Based on the comments and after a lot of trial and error plus reading carefully the Miepython documentation (!), I managed to answer my question:</p>
<p>For the problem that I wasn't able to reproduce scattering and backscatter efficiencies <span class="math-container">$qsca$</span> and <span class="math-container">$qback$</span> (or corresponding cross sections <span class="math-container">$\sigma$</span> in the first plot from values of the phase function in the second plot (update) in my question: I had two problems, or better say there were two pitfalls:</p>
<ul>
<li>I made an equidistant <span class="math-container">$\vartheta$</span>-vector and converted to <span class="math-container">$mu = cos(\vartheta)$</span>, calculated the phase function <span class="math-container">$p(\mu)$</span> and plotted it as polar plot and histogram, which I tried then to integrate (sum up). That was of course not working as <span class="math-container">$\mu = cos(\vartheta)$</span> is no longer equidistant. Now I'm doing an equidistant <span class="math-container">$\mu=cos(\vartheta)$</span>-vector for calculation and integration. For plotting against <span class="math-container">$\vartheta$</span> I just calculate <span class="math-container">$\vartheta$</span> from <span class="math-container">$\mu$</span> using <span class="math-container">$arccos()$</span>. The result is non-equidistant datapoints in the histogram then, but when the resolution is high enough, you don't see it..</li>
<li>There is a normalisation of the phase function in miepython, which is by default set to the single scattering albedo <span class="math-container">$a = qsca/qext$</span>. So after summing up (integrating), to get <span class="math-container">$qsca$</span> you need to multiply by <span class="math-container">$qext$</span>. (But there are other options you can pass to Miepython for normalisation, the default is just <span class="math-container">$qsca/qext$</span>... see <a href="https://miepython.readthedocs.io/en/latest/03a_normalization.html" rel="nofollow noreferrer">https://miepython.readthedocs.io/en/latest/03a_normalization.html</a>)</li>
</ul>
<p>With this, I could reproduce values for <span class="math-container">$qsca$</span> and <span class="math-container">$qback$</span> at 500 nm. Here is the code:</p>
<pre><code>import numpy as np
import miepython
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({'font.size': 10})

###############################################################################
# Parameters
###############################################################################
m = 1.335 - 1E-9j  # complex refractive index for water droplet at 500 nm from https://refractiveindex.info/?shelf=main&amp;book=H2O&amp;page=Daimon - 24.0C
lambda0 = 0.5 # in mu, i.e. 500 nm
diameter = 1 # in mu
r = 0.5*diameter  # radius of aerosol
x = 2*np.pi*r/lambda0 # size parameter

###############################################################################
# Miepython calculation
# - Calculate scattering efficiencies
# - Calculate phase function and plot as polar plot and histogram
###############################################################################
mu = np.linspace(1,-1,100001) # make mu eqidistant and calculate (for the plots) theta  from it, not vice versa (important as we integrate later the histogram, which is a function of mu, so mu needs to be equidistant for easy sum up -&gt; integrating)
theta = np.arccos(mu)*180/np.pi

# Calculate extinction, scattering and backscatter efficiency
qext500, qsca500, qback500, g2_500 = miepython.mie(m,x)

# Now calculate the phase function
normalisation='albedo' # default is albedo a = qsca/qext, see https://miepython.readthedocs.io/en/latest/03a_normalization.html
p_mu = miepython.i_unpolarized(m, x, mu, normalisation) # p_mu is now normalised to a = qsca/qext
# Test alternative way of calculating the phase function:
s1, s2 = miepython.mie_S1_S2(m, x, mu)
phase = (np.abs(s1)**2 + np.abs(s2)**2 )/2
if np.allclose(p_mu, phase):
    print('p_mu and phase are the same vectors')

# plot phase function as polar plot and histogram
fig, (ax1, ax2) = plt.subplots(1,2,figsize=(12,5))
ax1=plt.subplot(121, projection='polar')
ax1.plot(theta/180*np.pi, p_mu, color='blue')
ax1.set_title(r'm = ' + str(np.around(m,9)) + ', d = '+str(diameter) + r' <span class="math-container">$\mu$</span>m, <span class="math-container">$\lambda$</span> = ' + str(np.around(1000*lambda0,0)) + r' nm')

ax2 = plt.subplot(122)
ax2.plot(theta, p_mu)
ax2.set_xlabel(r'<span class="math-container">$\vartheta$</span> (<span class="math-container">$^\circ$</span>)')
ax2.set_ylabel('Unpolarized Scattered light [1/sr]')
ax2.grid(which = 'minor', color='grey', linestyle=':')
ax2.grid(which = 'major', color='grey', linestyle=':')
ax2.set_title(r'm = ' + str(np.around(m,9)) + ', d = '+str(diameter) + r' <span class="math-container">$\mu$</span>m, <span class="math-container">$\lambda$</span> = ' + str(np.around(1000*lambda0,0)) + r' nm')

plt.tight_layout()
plt.show()

fig.savefig('phase function.png', dpi=fig.dpi, bbox_inches='tight')
  
###############################################################################
# Now integrate histogram to reproduce scattering and backscatter effiencies 
###############################################################################
# print values for comparison
print('qback500 = ' + str(qback500))
print('qsca500 = ' + str(qsca500))
print('qext500 = ' + str(qext500))
print('qsca and qext is almost the same, because of almost no absorption (imaginary part of complex refractive index, so extinction (loss) is mainly caused by scattering...)')

print('\nNow calculating the scattering efficency from the phase function:')

# spacing dmu for integration
oben = np.abs(mu[1:])
unten = np.abs(mu[:-1])
dmu = np.mean(np.abs(oben - unten))
print('\ndmu = ' + str(dmu))
print('dmu is the spacing on the mu = cos(theta) axis')

# qsca from histogram:
print('Integration of phase function over full solid angle:')
print('- There is no dependency on the azimuth, so integration over azimuth gives 2*pi and the remaining integration is over mu only.')
print('- The integration over mu is from 1 to -1, i.e. from theta=0° to 180°')
integral_p_mu = 2*np.pi*dmu*np.sum(p_mu)
print('Integral over p_mu histogram = ' + str(integral_p_mu))
# To get qsca from this, remind that p_mu was normalised to a = qsca/qext. Thus, to get qsca, one needs to multiply the obtained result here with qext:
print('\nqsca calculated from phase function = ' + str(integral_p_mu * qext500))

print('\nBackscattering efficency calculated from physe function')
idx180 = np.where(np.abs(180 - theta) &lt; 0.01)[0][0] # index of backward direction
# Again, we need to multiply by qext because of the normalisation 
# We use only the 180° datapoint of the phase function, but still need to integrate over the solid angle, which just gives 4*pi')
print('p_mu[idx180]*qext500*4*np.pi = ' + str(p_mu[idx180]*qext500*4*np.pi))
</code></pre>
<p>And the answer to my original question is:</p>
<ul>
<li><span class="math-container">$qsca$</span> is the total scattering efficiency in all directions (integration over a full sphere of <span class="math-container">$4\pi$</span>), i.e. not only forward half-dome or so..</li>
<li><span class="math-container">$qback$</span> is scattering in the backward direction of <span class="math-container">$180^\circ$</span> only (one single direction, not the entire rear half-dome or so...)</li>
</ul>

CC BY-SA 4.0
2023-03-30T16:02:52.003Z
Dear @gansub Yes, many thanks for your help!
Charles Tucker 3
CC BY-SA 4.0
2023-03-30T07:25:54.307Z
2023-03-30T23:51:45.483Z
2023-03-30T23:51:45.483Z
Charles Tucker 3
<p>I’m usually on Physics SE, but this question seems to fit better into Earth Sciences:</p>
<p>I’m interested in radiation transfer and in particular Mie scattering in the atmosphere. So I did some calculation with Miepython (<a href="https://pypi.org/project/miepython/" rel="nofollow noreferrer">https://pypi.org/project/miepython/</a>). Below is an example result for extinction, scattering and backscattering on small water droplets of 1000 nm diameter.</p>
<p><a href="https://i.sstatic.net/JtSsv.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/JtSsv.png" alt="enter image description here" /></a></p>
<p>Unfortunately, I’m a bit puzzled about the results.</p>
<p>First, is there any common definition of the angular range in which scattering and backscattering is? The Miepython manual doesn’t state anything about. E.g., is scattering maybe in 0<span class="math-container">$^\circ$</span>-180<span class="math-container">$^\circ$</span> and backscattering in 90<span class="math-container">$^\circ$</span>-180<span class="math-container">$^\circ$</span>?</p>
<p>Second, the backscattering shows a lot of variations with wavelength. If I increase the particle size, these get even more. What is the physical explanation for this zig-zag behaviour of alternating maxima and minima?</p>
<p><strong>Update</strong></p>
<p>After receiving helpful comments I calculated a polar plot, following instructions from <a href="https://miepython.readthedocs.io/en/latest/03_angular_scattering.html" rel="nofollow noreferrer">https://miepython.readthedocs.io/en/latest/03_angular_scattering.html</a> (incl. axes labeling...) for unpolarized radiation.</p>
<p>The result is:</p>
<p><a href="https://i.sstatic.net/9DAql.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/9DAql.png" alt="enter image description here" /></a></p>
<p>However, I tried to integrate (sum up) over different ranges and multiplied by the geometric cross section (what miepython requires you to do, and which has been done already for the plots above), but couldn't reproduce numbers of the above plot at 500 nm.</p>

CC BY-SA 4.0
2022-07-07T12:16:52.373Z
@gansub Very good suggestion, I updated the question now with a polar plot. For the maxima and mimima question: Do you have a paper or some link for further reading explaining this in more detail?
Charles Tucker 3
CC BY-SA 4.0
2022-07-07T15:51:37.767Z
@gansub Alright, but the unit on the polar plot seems to be 1/sr. In which direction is the cone? Suppose I have a lidar and know its angular field of view. Do I just read the value for 180 deg from the polar plot and multiply with my fov to get the backscatter? Or do I need to calculate the solid angle of my aperture from the distance of the scattering height?
Charles Tucker 3
CC BY-SA 4.0
2022-07-10T14:40:40.410Z
@gansub I got the book now, but page 385 seems to be wrong in my version. Could you point me to the section number? For the other question: How do I get from the polar plot of the phase function to the backscattering or scattering cross-section (my first plot)? Thought I just need to integrate and divide by the geometric cross section, but that didn't reproduce numbers in the first plot (tested exemplarily for lambda = 500 nm)...
Charles Tucker 3
<atmospheric-radiation><radiative-transfer>
Interpretation of Mie scattering
Charles Tucker 3
# Interpretation of Mie scattering

(Asked by: username_0 on 2022-07-07)

I’m usually on Physics SE, but this question seems to fit better into Earth Sciences:

I’m interested in radiation transfer and in particular Mie scattering in the atmosphere. So I did some calculation with Miepython (<https://pypi.org/project/miepython/>). Below is an example result for extinction, scattering and backscattering on small water droplets of 1000 nm diameter.

[![enter image description here](https://i.sstatic.net/JtSsv.png)](https://i.sstatic.net/JtSsv.png)

Unfortunately, I’m a bit puzzled about the results.

First, is there any common definition of the angular range in which scattering and backscattering is? The Miepython manual doesn’t state anything about. E.g., is scattering maybe in 0$^\circ$-180$^\circ$ and backscattering in 90$^\circ$-180$^\circ$?

Second, the backscattering shows a lot of variations with wavelength. If I increase the particle size, these get even more. What is the physical explanation for this zig-zag behaviour of alternating maxima and minima?

**Update**

After receiving helpful comments I calculated a polar plot, following instructions from <https://miepython.readthedocs.io/en/latest/03_angular_scattering.html> (incl. axes labeling...) for unpolarized radiation.

The result is:

[![enter image description here](https://i.sstatic.net/9DAql.png)](https://i.sstatic.net/9DAql.png)

However, I tried to integrate (sum up) over different ranges and multiplied by the geometric cross section (what miepython requires you to do, and which has been done already for the plots above), but couldn't reproduce numbers of the above plot at 500 nm.

**Tags:** atmospheric-radiation, radiative-transfer

> **Comments:**
> **username_0**: @gansub Very good suggestion, I updated the question now with a polar plot. For the maxima and mimima question: Do you have a paper or some link for further reading explaining this in more detail?  
> **username_0**: @gansub Alright, but the unit on the polar plot seems to be 1/sr. In which direction is the cone? Suppose I have a lidar and know its angular field of view. Do I just read the value for 180 deg from the polar plot and multiply with my fov to get the backscatter? Or do I need to calculate the solid angle of my aperture from the distance of the scattering height?  
> **username_0**: @gansub I got the book now, but page 385 seems to be wrong in my version. Could you point me to the section number? For the other question: How do I get from the polar plot of the phase function to the backscattering or scattering cross-section (my first plot)? Thought I just need to integrate and divide by the geometric cross section, but that didn't reproduce numbers in the first plot (tested exemplarily for lambda = 500 nm)...  

---

## 1 answers

### Answer by username_0 on 2023-03-30. Score: 4

Based on the comments and after a lot of trial and error plus reading carefully the Miepython documentation (!), I managed to answer my question:

For the problem that I wasn't able to reproduce scattering and backscatter efficiencies $qsca$ and $qback$ (or corresponding cross sections $\sigma$ in the first plot from values of the phase function in the second plot (update) in my question: I had two problems, or better say there were two pitfalls:

* I made an equidistant $\vartheta$-vector and converted to $mu = cos(\vartheta)$, calculated the phase function $p(\mu)$ and plotted it as polar plot and histogram, which I tried then to integrate (sum up). That was of course not working as $\mu = cos(\vartheta)$ is no longer equidistant. Now I'm doing an equidistant $\mu=cos(\vartheta)$-vector for calculation and integration. For plotting against $\vartheta$ I just calculate $\vartheta$ from $\mu$ using $arccos()$. The result is non-equidistant datapoints in the histogram then, but when the resolution is high enough, you don't see it..
* There is a normalisation of the phase function in miepython, which is by default set to the single scattering albedo $a = qsca/qext$. So after summing up (integrating), to get $qsca$ you need to multiply by $qext$. (But there are other options you can pass to Miepython for normalisation, the default is just $qsca/qext$... see <https://miepython.readthedocs.io/en/latest/03a_normalization.html>)

With this, I could reproduce values for $qsca$ and $qback$ at 500 nm. Here is the code:

```
import numpy as np
import miepython
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({'font.size': 10})

###############################################################################
# Parameters
###############################################################################
m = 1.335 - 1E-9j  # complex refractive index for water droplet at 500 nm from https://refractiveindex.info/?shelf=main&book=H2O&page=Daimon - 24.0C
lambda0 = 0.5 # in mu, i.e. 500 nm
diameter = 1 # in mu
r = 0.5*diameter  # radius of aerosol
x = 2*np.pi*r/lambda0 # size parameter

###############################################################################
# Miepython calculation
# - Calculate scattering efficiencies
# - Calculate phase function and plot as polar plot and histogram
###############################################################################
mu = np.linspace(1,-1,100001) # make mu eqidistant and calculate (for the plots) theta  from it, not vice versa (important as we integrate later the histogram, which is a function of mu, so mu needs to be equidistant for easy sum up -> integrating)
theta = np.arccos(mu)*180/np.pi

# Calculate extinction, scattering and backscatter efficiency
qext500, qsca500, qback500, g2_500 = miepython.mie(m,x)

# Now calculate the phase function
normalisation='albedo' # default is albedo a = qsca/qext, see https://miepython.readthedocs.io/en/latest/03a_normalization.html
p_mu = miepython.i_unpolarized(m, x, mu, normalisation) # p_mu is now normalised to a = qsca/qext
# Test alternative way of calculating the phase function:
s1, s2 = miepython.mie_S1_S2(m, x, mu)
phase = (np.abs(s1)**2 + np.abs(s2)**2 )/2
if np.allclose(p_mu, phase):
    print('p_mu and phase are the same vectors')

# plot phase function as polar plot and histogram
fig, (ax1, ax2) = plt.subplots(1,2,figsize=(12,5))
ax1=plt.subplot(121, projection='polar')
ax1.plot(theta/180*np.pi, p_mu, color='blue')
ax1.set_title(r'm = ' + str(np.around(m,9)) + ', d = '+str(diameter) + r' $\mu$m, $\lambda$ = ' + str(np.around(1000*lambda0,0)) + r' nm')

ax2 = plt.subplot(122)
ax2.plot(theta, p_mu)
ax2.set_xlabel(r'$\vartheta$ ($^\circ$)')
ax2.set_ylabel('Unpolarized Scattered light [1/sr]')
ax2.grid(which = 'minor', color='grey', linestyle=':')
ax2.grid(which = 'major', color='grey', linestyle=':')
ax2.set_title(r'm = ' + str(np.around(m,9)) + ', d = '+str(diameter) + r' $\mu$m, $\lambda$ = ' + str(np.around(1000*lambda0,0)) + r' nm')

plt.tight_layout()
plt.show()

fig.savefig('phase function.png', dpi=fig.dpi, bbox_inches='tight')
  
###############################################################################
# Now integrate histogram to reproduce scattering and backscatter effiencies 
###############################################################################
# print values for comparison
print('qback500 = ' + str(qback500))
print('qsca500 = ' + str(qsca500))
print('qext500 = ' + str(qext500))
print('qsca and qext is almost the same, because of almost no absorption (imaginary part of complex refractive index, so extinction (loss) is mainly caused by scattering...)')

print('\nNow calculating the scattering efficency from the phase function:')

# spacing dmu for integration
oben = np.abs(mu[1:])
unten = np.abs(mu[:-1])
dmu = np.mean(np.abs(oben - unten))
print('\ndmu = ' + str(dmu))
print('dmu is the spacing on the mu = cos(theta) axis')

# qsca from histogram:
print('Integration of phase function over full solid angle:')
print('- There is no dependency on the azimuth, so integration over azimuth gives 2*pi and the remaining integration is over mu only.')
print('- The integration over mu is from 1 to -1, i.e. from theta=0° to 180°')
integral_p_mu = 2*np.pi*dmu*np.sum(p_mu)
print('Integral over p_mu histogram = ' + str(integral_p_mu))
# To get qsca from this, remind that p_mu was normalised to a = qsca/qext. Thus, to get qsca, one needs to multiply the obtained result here with qext:
print('\nqsca calculated from phase function = ' + str(integral_p_mu * qext500))

print('\nBackscattering efficency calculated from physe function')
idx180 = np.where(np.abs(180 - theta) < 0.01)[0][0] # index of backward direction
# Again, we need to multiply by qext because of the normalisation 
# We use only the 180° datapoint of the phase function, but still need to integrate over the solid angle, which just gives 4*pi')
print('p_mu[idx180]*qext500*4*np.pi = ' + str(p_mu[idx180]*qext500*4*np.pi))

```

And the answer to my original question is:

* $qsca$ is the total scattering efficiency in all directions (integration over a full sphere of $4\pi$), i.e. not only forward half-dome or so..
* $qback$ is scattering in the backward direction of $180^\circ$ only (one single direction, not the entire rear half-dome or so...)

> **Comments:**
> **username_0**: Dear @gansub Yes, many thanks for your help!  

---

