coefficents calculation at large size parameter gives unexpected results

[ysBach]

The following code gives results which I didn't expect:

import numpy as np
import matplotlib.pyplot as plt

# if miepython is missing, do `pip install miepython`
import miepython as mp

m = 1.5 - 0.01j
x = np.logspace(-1, 5, 20)  # also in microns

qext, qsca, qback, g = mp.mie(m,x)
plt.plot(x, qext, 'b+-', label="Q_ext")
plt.plot(x, qsca, 'r+:', label="Q_sca")
plt.plot(x, qback, 'g+:', label="Q_back")
plt.plot(x, g, 'r+-', label="<cos>")
plt.plot(x, qext - g * qsca, 'kx-', label="Q_pr")
plt.xscale('log')
plt.yscale('log')
plt.xlabel("size parameter")
plt.legend()
plt.grid()

The resulting figure:

At size parameter larger than ~ 5000, all coefficients except for qext (expected to be ~ 2) deviate largely, and qpr becomes negative.

As I am not so familiar with scattering theory, I couldn't understand/interpret this result. Is this an expected behavior?

[scottprahl]

You definitely found a bug! Nice work.

My excuse is that I have not cared about very large spheres and therefore just did not test these cases. I will upload a new jupyter notebook to the doc section showing

It turns out that the flag to determine if upwards or downward recurrence relations should be used was messed up. I have a fix for this that solves the Qsca, Qext, and Qpr problems. I am still a bit unsure if the fixes the problem with Qback, but I don't have the time to investigate it at the moment.

I uploaded the fix to miepython.py to github. Once I convince myself that Qback is also fine for large spheres, then I will release a new version.

Here is what your code generates with the new fix

[ysBach]

@scottprahl Thank you so much for your prompt reply and update of the package! :)

I have tested 273f46ff53ef4f00ca68d0723c87fa33100a2a8f and it works perfectly, at least seemingly.

FYI) The reason I used your miepython to my research is I wanted to calculate the black-body spectrum averaged $ Q\mathrm{pr} $, $ \bar{Q} $, of certain mineral, and now the results are as I expected: (The codes are very long and time consuming so I didn't attach any here) Using the previous code, the $\bar{Q}\mathrm{pr}$ value gets negative when $a > 100 \mathrm{\mu m}$, and the trend does not seem to converge to 1 even when particle size is very large.

But this new code shows that $\bar{Q}_\mathrm{pr}$ is always positive and converges to unity when the particle is very large (compared to the wavelengths of interest).

(Update: Furthermore, the shape of each curve and the trend of its change as a function of black-body temperature are also similar to other minerals, which is as I expected.)

[scottprahl]

Lovely graphs!

I want to thank you again for your bug report. It was just the right thing that I needed to see to track down the bug. I knew something was awry with large sphere calculations but attributed the problems to round-off when summing many terms. Your first graph showed me that the transition was sudden and therefore I was able to change one line of code that fixed calculations for all large sphere.