First point calculated by fieldnlay is always inf&nan

[10110111]

The following command always fails to calculate field for the first point in the output, at least when |z|≥300; both in multiprecision and double precision versions. For double precision z can be as low as 100 to reproduce the problem.

$ ./fieldnlay-mp -l 1 113 1.3346 0 -p -50 50 2 -50 50 2 300 300 1
Econv:0.15942 Hconv:0.170959
Econv:0.120661 Hconv:0.128241
Econv:0.0964606 Hconv:0.0910241
Field evaluation failed to converge an nmax = 189
Near-field convergence threshold: 1e-14
         X,          Y,          Z,         Ex.r,         Ex.i,         Ey.r,         Ey.i,         Ez.r,         Ez.i,         Hx.r,         Hx.i,         Hy.r,         Hy.i,         Hz.r,         Hz.i
-50.0000000, -50.0000000, 300.0000000, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan
-50.0000000, 50.0000000, 300.0000000, -3.69341e-01, -4.20365e-01, +1.84654e-02, -7.65961e-03, -4.59793e-02, -2.67651e-01, +4.54062e-05, -1.23898e-04, -9.81315e-04, -1.12519e-03, +1.17023e-04, +6.59271e-04
50.0000000, -50.0000000, 300.0000000, -3.69341e-01, -4.20365e-01, +1.84654e-02, -7.65961e-03, +4.59793e-02, +2.67651e-01, +4.54062e-05, -1.23898e-04, -9.81315e-04, -1.12519e-03, -1.17023e-04, -6.59271e-04
50.0000000, 50.0000000, 300.0000000, -3.69341e-01, -4.20365e-01, -1.84654e-02, +7.65961e-03, +4.59793e-02, +2.67651e-01, -4.54062e-05, +1.23898e-04, -9.81315e-04, -1.12519e-03, +1.17023e-04, +6.59271e-04

Repeating calculation at the same point yields a finite value on the second attempt:

$ ./fieldnlay-mp -l 1 113 1.3346 0 -p -50 50 2 -50 50 2 300 300 2
Econv:0.15942 Hconv:0.170959
Econv:0.120661 Hconv:0.128241
Econv:0.0964606 Hconv:0.0910241
Field evaluation failed to converge an nmax = 189
Near-field convergence threshold: 1e-14
         X,          Y,          Z,         Ex.r,         Ex.i,         Ey.r,         Ey.i,         Ez.r,         Ez.i,         Hx.r,         Hx.i,         Hy.r,         Hy.i,         Hz.r,         Hz.i
-50.0000000, -50.0000000, 300.0000000, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan
-50.0000000, -50.0000000, 300.0000000, -3.69341e-01, -4.20365e-01, -1.84654e-02, +7.65961e-03, -4.59793e-02, -2.67651e-01, -4.54062e-05, +1.23898e-04, -9.81315e-04, -1.12519e-03, -1.17023e-04, -6.59271e-04
-50.0000000, 50.0000000, 300.0000000, -3.69341e-01, -4.20365e-01, +1.84654e-02, -7.65961e-03, -4.59793e-02, -2.67651e-01, +4.54062e-05, -1.23898e-04, -9.81315e-04, -1.12519e-03, +1.17023e-04, +6.59271e-04
-50.0000000, 50.0000000, 300.0000000, -3.69341e-01, -4.20365e-01, +1.84654e-02, -7.65961e-03, -4.59793e-02, -2.67651e-01, +4.54062e-05, -1.23898e-04, -9.81315e-04, -1.12519e-03, +1.17023e-04, +6.59271e-04
50.0000000, -50.0000000, 300.0000000, -3.69341e-01, -4.20365e-01, +1.84654e-02, -7.65961e-03, +4.59793e-02, +2.67651e-01, +4.54062e-05, -1.23898e-04, -9.81315e-04, -1.12519e-03, -1.17023e-04, -6.59271e-04
50.0000000, -50.0000000, 300.0000000, -3.69341e-01, -4.20365e-01, +1.84654e-02, -7.65961e-03, +4.59793e-02, +2.67651e-01, +4.54062e-05, -1.23898e-04, -9.81315e-04, -1.12519e-03, -1.17023e-04, -6.59271e-04
50.0000000, 50.0000000, 300.0000000, -3.69341e-01, -4.20365e-01, -1.84654e-02, +7.65961e-03, +4.59793e-02, +2.67651e-01, -4.54062e-05, +1.23898e-04, -9.81315e-04, -1.12519e-03, +1.17023e-04, +6.59271e-04
50.0000000, 50.0000000, 300.0000000, -3.69341e-01, -4.20365e-01, -1.84654e-02, +7.65961e-03, +4.59793e-02, +2.67651e-01, -4.54062e-05, +1.23898e-04, -9.81315e-04, -1.12519e-03, +1.17023e-04, +6.59271e-04

Tested on a 32-bit x86 Linux system, with glibc 2.27 and GCC 7.5.0.

[kostyfisik]

Thank you for reporting! It looks like a real bug in the current master branch. You can check out to tagged release v2.2.0, and it provides reproducible results. However, there was no explicit convergence check, we were blindly using a Wiscombe criteria, which is now known not to be valid for near-field computations. So, while it will give you some value, most likely it is wrong for large values of size-parameters\distance from sphere\refractive index (especially imaginary part). You might be interested in related closed issue #19

$ ./fieldnlay-mp -l 1 113 1.3346 0 -p -50 50 2 -50 50 2 300 300 2
         X,          Y,          Z,         Ex.r,         Ex.i,         Ey.r,         Ey.i,         Ez.r,         Ez.i,         Hx.r,         Hx.i,         Hy.r,         Hy.i,         Hz.r,         Hz.i
-50.0000000, -50.0000000, 300.0000000, +6.01901e-01, +3.17665e-01, +1.63639e-01, -1.52411e-01, +6.84459e-02, +9.98730e-02, +4.18473e-04, -4.32869e-04, +1.91833e-03, +9.47548e-05, +2.81127e-04, +5.25052e-05
-50.0000000, -50.0000000, 300.0000000, +6.01901e-01, +3.17665e-01, +1.63639e-01, -1.52411e-01, +6.84459e-02, +9.98730e-02, +4.18473e-04, -4.32869e-04, +1.91833e-03, +9.47548e-05, +2.81127e-04, +5.25052e-05
-50.0000000, 50.0000000, 300.0000000, +6.01901e-01, +3.17665e-01, -1.63639e-01, +1.52411e-01, +6.84459e-02, +9.98730e-02, -4.18473e-04, +4.32869e-04, +1.91833e-03, +9.47548e-05, -2.81127e-04, -5.25052e-05
-50.0000000, 50.0000000, 300.0000000, +6.01901e-01, +3.17665e-01, -1.63639e-01, +1.52411e-01, +6.84459e-02, +9.98730e-02, -4.18473e-04, +4.32869e-04, +1.91833e-03, +9.47548e-05, -2.81127e-04, -5.25052e-05
50.0000000, -50.0000000, 300.0000000, +6.01901e-01, +3.17665e-01, -1.63639e-01, +1.52411e-01, -6.84459e-02, -9.98730e-02, -4.18473e-04, +4.32869e-04, +1.91833e-03, +9.47548e-05, +2.81127e-04, +5.25052e-05
50.0000000, -50.0000000, 300.0000000, +6.01901e-01, +3.17665e-01, -1.63639e-01, +1.52411e-01, -6.84459e-02, -9.98730e-02, -4.18473e-04, +4.32869e-04, +1.91833e-03, +9.47548e-05, +2.81127e-04, +5.25052e-05
50.0000000, 50.0000000, 300.0000000, +6.01901e-01, +3.17665e-01, +1.63639e-01, -1.52411e-01, -6.84459e-02, -9.98730e-02, +4.18473e-04, -4.32869e-04, +1.91833e-03, +9.47548e-05, -2.81127e-04, -5.25052e-05
50.0000000, 50.0000000, 300.0000000, +6.01901e-01, +3.17665e-01, +1.63639e-01, -1.52411e-01, -6.84459e-02, -9.98730e-02, +4.18473e-04, -4.32869e-04, +1.91833e-03, +9.47548e-05, -2.81127e-04, -5.25052e-05

[kostyfisik]

I've updated the convergence check, it is available in quasar branch now. The early convergence happens if the change of all components of E or H fields is <= of previously evaluated E or H component value multiplied by 1e-14 sequentially for two terms of Mie series. This helps to avoid false convergence when accidentally taking near-field point at some zero of vector spherical harmonic. I have also added incident field evaluation as it was proposed in #44, thank you very much for this proposal! The updated result:

$ ./fieldnlay-dp -l 1 113 1.3346 0 -p -50 50 2 -50 50 2 300 300 1
Near-field early convergence at nmax = 146
Near-field early convergence at nmax = 146
Near-field early convergence at nmax = 146
Near-field early convergence at nmax = 146
Number of multipoles used in Mie series nmax=183
         X,          Y,          Z,         Ex.r,         Ex.i,         Ey.r,         Ey.i,         Ez.r,         Ez.i,         Hx.r,         Hx.i,         Hy.r,         Hy.i,         Hz.r,         Hz.i
-50.0000000, -50.0000000, 300.0000000, -3.59081e-01, -5.23959e-01, +1.55091e-02, +2.24917e-02, -8.22586e-02, -2.42369e-01, +4.47763e-05, +1.63268e-04, -9.54081e-04, -1.40017e-03, -2.13323e-04, -5.92162e-04
-50.0000000, 50.0000000, 300.0000000, -3.59081e-01, -5.23959e-01, -1.55091e-02, -2.24917e-02, -8.22586e-02, -2.42369e-01, -4.47763e-05, -1.63268e-04, -9.54081e-04, -1.40017e-03, +2.13323e-04, +5.92162e-04
50.0000000, -50.0000000, 300.0000000, -3.59081e-01, -5.23959e-01, -1.55091e-02, -2.24917e-02, +8.22586e-02, +2.42369e-01, -4.47763e-05, -1.63268e-04, -9.54081e-04, -1.40017e-03, -2.13323e-04, -5.92162e-04
50.0000000, 50.0000000, 300.0000000, -3.59081e-01, -5.23959e-01, +1.55091e-02, +2.24917e-02, +8.22586e-02, +2.42369e-01, +4.47763e-05, +1.63268e-04, -9.54081e-04, -1.40017e-03, +2.13323e-04, +5.92162e-04

See more details here https://github.com/ovidiopr/scattnlay/blob/quasar/src/nmie-nearfield.hpp#L322-L329

[kostyfisik]

@10110111 BTW, I would like to warn you once more. Typical usage of Scattnaly is for nanophotonics and plasmonics, where the most important are just few multipoles in Mie series. It is rather rear to see nmax> 30. So, in your case with nmax>100 is far away from typical usage.

Choosing correctly cutoff value of nmax is a hard problem, which is probably not yet reliably solved yet. Presently I set the maximum value choosing between Wiscombe limit, Le Ru limit [1], and my own assumption (the last one is mostly for the field outside the sphere, and it depends on near-field evaluation point distance from the sphere and based on Le Ru). Additionally, I have recently added extended evaluation of special functions of maximum order using Kapteyn inequality [2]. For near field outside the particle, it should be rather stable now. For the field inside the bulk particle it should be rather stable too, for multilayered particles we plan to rewrite the equations in a more numerically stable form (this part is not yet implemented).
[1] J. R. Allardice and E. C. Le Ru, Convergence of Mie Theory Series: Criteria for Far-Field and near-Field Properties, Appl. Opt. 53, 7224 (2014).
[2]
H. Du, Mie-Scattering Calculation, Appl. Opt. 43, 1951 (2004).

I would also recommend you to try Ilia Rasskazov code Stratify https://gitlab.com/iliarasskazov/stratify It is a Matlab implementation, I used it to cross-check our near-field implementation and it was good. Another possible option is a code by Honoh Suzuki [3], which has a multiprecision version. At least, I am planning to add a lot of tests of near-field evaluation using these two codes later on.
[3]
H. Suzuki and I.-Y. Sandy Lee, Mie Scattering Field inside and near a Coated Sphere: Computation and Biomedical Applications, Journal of Quantitative Spectroscopy and Radiative Transfer 126, 56 (2013).

You might also be interested in http://www.philiplaven.com/mieplot.htm code, which is mostly written with your case in mind - optically large water droplets..

[10110111]

Thanks for the update, looks great!

It is rather rear to see nmax> 30. So, in your case with nmax>100 is far away from typical usage.

But non-typical doesn't mean unsupported, does it?

As for the recommendations, from what I have currently tried for near-field calculations, Scattnlay appears to be the most stable and easy to use package, without dependency on paid/closed-source software like MATLAB.

The Suzuki paper doesn't seem to give any link to the actual code, so I don't know whether it's any good for me.

MiePlot seems to only handle far field, and it's closed-source. For far field, aside from Scattnlay, there's miepython which works even for large spheres, but doesn't support multiprecision. So even here Scattnlay seems to be the best :)

[kostyfisik]

In general, I am interested to make our code to be universal one, limited only with available computational resources. For scattering it seems to work OK now, I have tested basic functions to be evaluated correctly for nmax of several thousands. I have also tested large water droplets available from the papers (still needs to add this data as test-cases). However, near-field is way less explored. Could you create a separate issue and to describe in more details your setup? E.g. what is the range for size of your droplets? Which part of spectrum you need? Are you interested in bulk or core-shell (e.g. soot aerosol)? What are the refractive indices of materials (essential, how much absorption do they have), etc. ?

Suzuki code is available at https://scattport.org/index.php/light-scattering-software/mie-type-codes/list/521-bhfield, I also attach it here. BTW, another good option can be https://github.com/gevero/py_gmm (I never used it, but Giovanni is a real expert in this field, so I expect his code to be good (as far as I understand, it works well for large particles without mulitprecision).

bhfield-121005.zip

[10110111]

Could you create a separate issue and to describe in more details your setup?

Well, I could, but I'm not sure what that issue would be about. Aside from #44, #43 and this issue, I've not seen any problems with the near field calculations.

E.g. what is the range for size of your droplets? Which part of spectrum you need? Are you interested in bulk or core-shell (e.g. soot aerosol)? What are the refractive indices of materials (essential, how much absorption do they have), etc. ?

I'm considering 10 μm water drops (uncoated spheres with n≈1.33, no absorption) in a vacuum, illuminated by a pulse of light containing several semiperiods of oscillation, with wavelength of 650 nm. This is done via Fourier transform, which spans wavelengths from about 100 nm to about 1.3 mm (for a single semiperiod), which implies size parameters from 0.05 to 600.

This is a similar simulation to the one here, but for near field (spanning about 10 sphere radii) instead of the far field.

[kostyfisik]

So, it seems that all your issues are resolved now in quasar branch, I will check them once again after merging to master and will finally close them after that. To verify the correctness of the results for your case, I will need to create a mpmath Python version for near-field evaluation and create an integration test from your data. I already have the main part of the mpmath test generator (https://github.com/ovidiopr/scattnlay/blob/master/tests/mpmath_special_functions_test_generator.py), it works nicely for total scattering. Extending it to near-field and using your case as a target should not be too hard (and I was going to implement it anyway, I will simply extend my current set of test cases (https://github.com/ovidiopr/scattnlay/blob/master/tests/mpmath_input_arguments.py) with your example (note, that some of these tests are still failing in Scattnlay, so it is work in progress).

[kostyfisik]

Sorry, I've just checked that test for water droplet with size parameter 10 000 passes OK in double precision for scattering, so good chances that for near field there is no need in multiple precision for your case too. https://github.com/ovidiopr/scattnlay/blob/master/tests/test_bulk_sphere.cc#L52

[kostyfisik]

@10110111 FYI, I've put a 10 mkm water droplet to Mie calculator webapp (it is using the same C++ under the hood in quasar git branch), it seems to provide some reasonable results in area of interest. The app is available at https://physics.itmo.ru/en/mie#/nearfield , it took about 2 min to run the simulation on my laptop for the plot below (I've capped the field because of the hotspot behind the droplet).

[10110111]

Looks good. Nice Arago spot :)

I've capped the field because of the hotspot behind the droplet

It might be a good idea to use a logarithmic color scale, similar to the one used by Philip Laven in MiePlot, maybe with another cutoff value (image taken from here):