Bigspheremielens

[briandleahy]

This PR makes mielens more stable at large size parameters, by changing the default maximum number of terms in the Mie scattering solutions to include.

The old problem was that, for large enough spheres, the Mie scattering matrices were evaluating something of the form (0 - 0) / (0 - 0), which gave nans. To fix this, I found the maximum number of terms I could include in the series, as a function of size_parameter for a fixed refractive index ratio. Empirically, this grows ever so slightly faster than linearly. Next, I found the minimum number of terms I could include for a ~1e-6 fractional error, as a function of size_parameter. This grows linearly, with a sizeable gap for all size parameters between the minimum number of terms and the maximum, i.e. it should be always possible to compute this numerically, up to size_parameter=1e4 (or a 790 um sphere), which is the largest that I checked. Then I found a simple functional form for the number of terms in the series to include, as a function of the size parameter.

The old code used the first 4 * size_parameter terms in the series, which caused errors at sphere sizes of about 14 um. The new code uses the first 25 + 1.1 * size_parameter terms in the series, as determined above. Both of these are different than what Bohren and Huffman suggest, which is (I believe) size_parameter**(4/3). I used a different scaling because (i) my scaling works empirically, and (ii) using their scaling gives results that crash for large enough size_parameter.

All these empirical tests are local and obviously not part of holopy; I'm attaching the script and results below so you can take a look if you want. Be aware; it takes a while to run (hours on my machine, since computing Mie scattering matrices for large spheres is very slow).

I also only checked this with an index ratio of 1.2 (polystyrene in water); if we want, I can check this for a higher contrast. (Lower contrasts should be simpler.) fix-mielens-nans.tar.gz

[vnmanoharan]

How does the new scaling rule affect our calculation times for size parameters of order 1-10? At the lower end of this range it would seem that we're now calculating a lot more coefficients.

[briandleahy]

There is essentially no difference for size parameters of order 1-10.

Attached is a script which times results for just the mielensfunctions portion of calc_holo. Right now the actual computation of the hologram through mielensfunctions is only 50% of the time; the xarray calls outside, in calc_holo and ScatteringTheory, make up the other 50%.

If I time the results on the new branch, as

git checkout f583b687
python3 time_mielens.py

I get

Time per call, ka=1: 0.015 s Time per call, ka=10: 0.016 s

whereas if I run the results on the old branch

git checkout 1ab50bcf
python3 time_mieelns.py

I get

Time per call, ka=1: 0.012 s Time per call, ka=10: 0.016 s

This is because, for small size parameters, much more time in the calculation is spent doing the integrals than evaluating the scattering matrices.

Calculating a hologram through calc_holo adds about 20 ms of overhead time on my machine, which swamps any of the changes from small size parameters. time_mielens.tar.gz

[barkls]

Everything looks good! I do think it's worth checking that the scaling you chose works at a higher index contrast. You mentioned there's a wide gap between min/max viable number of terms, so there should be room to adjust it if necessary. That way there's no chance someone will have to repeat this analysis in the future!

[briandleahy]

I re-ran the analysis at a higher index contrast. There is essentially no difference between an index ratio of 1.2 (polystyrene in water) or an index ratio of 2.5 (n=3.32 in water).

If you want me to go even higher i can; just let me know.

[barkls]

Thanks for verifying that! I also appreciate the plot.