scottprahl / miepython

Mie scattering of light by perfect spheres
MIT License
163 stars 56 forks source link

Difference between your implementation and that of Maetzler #7

Closed Apollo3zehn closed 5 years ago

Apollo3zehn commented 5 years ago

When I use your online calculator or this python module with following settings:

wavelength = 1500e-9
r = 4.710E-07
k = 2*3.14/wavelength
x = k * r             # dimensionless Mie size parameter = 1.9719
m = 1.4100 - 0.0080j  # index of refraction of sphere

qext, qsca, qback, g = miepython.mie(m,x)

print("The qext efficiency  is %.3f" % qext)
print("The scattering efficiency  is %.3f" % qsca)
print("The backscatter efficiency is %.3f" % qback)
print("The scattering anisotropy  is %.3f" % g)

then I get the following results:

The qext efficiency  is 1.155
The scattering efficiency  is 1.090
The backscatter efficiency is 0.097
The scattering anisotropy  is 0.658

With the script of Maetzler (bottom of page https://omlc.org/software/mie/), if I do:

mie(complex(1.410E+00, -8.000E-03), 2*pi/1500e-9*4.710E-01*1e-6).'

I get:

ans =

    1.4100
   -0.0080
    1.9729
    1.0879
    1.1556
   -0.0677
    0.1161
    0.6533
    0.1005

where ans(4) = qext = 1.0879 and ans(5) = qsca = 1.1556.

So it seems that the values of the extinction coefficient and of the scattering coefficent are interchanged.

I digged deeper in your code and found that the Python Mie coefficients an and bn are different to that of Maetzler:

Maezler (an):

  0.310171419307873 - 0.473671739481655i
  0.037432108684527 - 0.200893621373526i
  0.000067151985416 - 0.021297737932288i
 -0.000019924584511 - 0.001285929209167i
 -0.000000825241069 - 0.000050906381267i
 -0.000000022452476 - 0.000001402834546i
 -0.000000000449367 - 0.000000028342714i
 -0.000000000006890 - 0.000000000437321i
 -0.000000000000083 - 0.000000000005320i

Python (an):

  0.317746891588117 + 0.454752059193760i
  0.045298858385873 + 0.197467459081684i
  0.000838661081755 + 0.021263464324966i
  0.000023228420817 + 0.001285707362372i
  0.000000830332704 + 0.000050900760559i
  0.000000022453522 + 0.000001402654704i
  0.000000000449301 + 0.000000028338490i
  0.000000000006889 + 0.000000000437247i
  0.000000000000083 + 0.000000000005319i

Most numbers are equal but the signs are different. Same is true for bn values.

I also tried this Java Applet (http://www.lightscattering.de/MieCalc/eindex.html), which is really hard to get it running and it provides confusing results. Sometimes it provides the same results as the Python script, sometimes it there are clearly wrong values. Therefore I also tried "bhmie Matlab" from this site: http://scatterlib.wikidot.com/mie

With this script, I get the same results as with the script of Maetzler.

This is very confusing. Do you know more about this issue?

Thanks Apollo3zehn

pisarik commented 5 years ago

Dear Apollo3zehn, I have run into the same problem. Eventually I switched to PyMieScatt, it gives consistent results with Maetzler code.

Apollo3zehn commented 5 years ago

Thanks for your link! There are much more Mie implementations than I was aware of a week ago. But hopefully in the end they all provide the same results :) At least with the same sign and only scaled by a constant factor.

scottprahl commented 5 years ago

Bohren and Huffman use positive values for the imaginary part of the complex index of refraction — see equation 2.48 in their book. Based on your observations, Maetzler seems to the same.

miepython follows the convention used by Wiscombe and assumes absorbing particles have a negative imaginary refractive index.

Try rerunning the Matlab code using positive values for the imaginary index of refraction. I expect that you'll get a much better match.

scottprahl commented 5 years ago

I forgot to mention, normalization for scattering functions varies considerably across codes. The miepython scattering function assumes that the integral of the scattering function over all 4pi steradians is equal to the single scattering albedo (see detailed documentation)

Apollo3zehn commented 5 years ago

Thank you, this solves two issues at once :) With inversed sign, the end results are now consistent (altough the sign of the complex part of an is still inversed, but this may be an implementation detail).

I found this publication to be very informative (for others who might stumble across this): Sign Conventions in Electromagnetic (EM) Waves.

scottprahl commented 5 years ago

Yes, the signs of the imaginary parts of A_n and B_n differ in code bases that make different assumptions about the imaginary part of the refraction. I have never seen any place that this affects and observable quantity. It is just an implementation detail :-)