Is it possible to use py_gmm for dielectric objects?

[paulmueller]

I would like to compute the near field behind a dielectric sphere using py_gmm. From the ipython notebooks, I can tell that I need to define a material in the epsilon folder (https://github.com/gevero/py_gmm/tree/master/epsilon).

Is there a simpler way of computing the field behind a sphere given just the radius in wavelengths and refractive indices of the sphere and the medium?

I appreciate any help to get me started.

Thanks, Paul

[gevero]

Hi Paul

You can certainly do that. Let me comment on a few things in order to help you out:

• even if the primary target of py_gmm are metals, it can work perfectly well with dielectric, as you can se here. Tabulated optical constants for metals, and the utilities to import them in py_gmm, are only provided for the user convenience.
• Ultimately, the optical constants are passed to py_gmm with the m_eps numpy array. This array has n_sphere rows, with n_sphere>2 and 2 columns, the first for the real part of the dielectric function, the second for the imaginary part of the dielectric function. All you have to do is to create the array, and fill the first column with the square of the refractive index of your sphere. You can leave the second column empty if your refractive index is purely real.
• All the other parameters are also simply passed as arrays, i.e. v_r contains the radii for all the spheres, m_xyz the (x,y,z) coordinates for all the spheres, and so on. If this is to complex, I would be happy to upload an example, provided that you give me some toy parameters for the systems you want to study.

Best

Giovanni

PS Please remember that py_gmm needs at least two spheres. If you want to study a single sphere you can simply add a very small dummy sphere very far away from the sphere of interest

[paulmueller]

Hi Giovanni,

thank you. It would very much help me if you could create a little toy example and if it is not too much trouble for you. I am actually looking for Mie scattering code that can also deal with large dielectric spheres. I have played around with GMM-FIELD, but there the size of the sphere is limited.

Is it, in principle, possible to use py_gmm with very large radii (e.g. 5-20 wavelengths) and a refractive index of glass (~1.46) embedded in a medium (~1.333)?

I am also looking into scattnlay (https://github.com/ovidiopr/scattnlay/issues/6#) and Kostya proposed to use BHFIELD (http://www.scattport.org/index.php/light-scattering-software/mie-type-codes/list/521-bhfield) for large spheres, because it has arbitrary accuracy.

Regards Paul

[gevero]

Hi Paul

Yes, in principle and in practice it is possible. I will happily put up a small toy example for you. It seems you need large silica spheres in water. In the meantime let me explain you how it works:

• First, remember that you need at least 2 spheres with py_gmm, so you can place your sphere of interest in the origin, and a much smaller sphere far away (let's say 100 radii away from the principal sphere)
• usually, a fully interacting calculation in the Generalized Multiparticle Mie framework requires much more memory that a simple Mie calculation. Morover the Wiscombe's criterion, which expresses the number of multipoles N to include in the calculations as a function of the size parameter x, (N~x+4.05x^(1/3)+2), does not apply
• Anyway, if you look at the AuAg example, you'll notice that among the many inputs these if an interaction parameter fint. If this parameter is set to 0, all the spheres are interacting, but if you set it to a number larger than (r1+r2)/d12, where r1 is the radius of the principal sphere, r2 is the radius of the dummy sphere, and d12 is the center to center distance, the spheres won't be interacting and the code will performe independent Mie calculations, with much lower memory requirements. In this context the code can handle sphere much larger than the wavelength, since I made sure that the computation of single sphere expansion coefficients is very stable, following this paper.

Finally, you should simply proceed as follows:

• Put your principal sphere in the origin
• Put your dummy sphere far away from the principal one and make it small, (radius 10 times smaller and 10-100 principal radii away from the main sphere)
• Set fint appropriately, so that the code performs simple mie calculation
• Chose the number of multipole expansion appropriately, in this case, since it is essentially a calculation for independent spheres, Wiscombe criterion should be ok.

best

Giovanni

[kostyfisik]

Dear Giovanni,

I am not sure, but I believe that the origin of Pauls problem is not physical (the number of multipoles needed) but a mathematical one (Bessels functions are oscillating). I have checked the downward recursive calculation of the Riccati-Bessel functions (it looks that in Scattnlay we evaluate them in a similar to GMM way) against Wolfram-Alpha, and it does not converge even for large number of elements in the recursion if the argument is large enough. (Note that Matlab`s functions give a little bit different result compared to Wolfram-Alpha, I believe that Matlab is wrong). So the obvious way is to increase the precision of the evaluation from 16 digits of double precision to any needed number of digits is to use some arbitrary precision lib, that why I mentioned BHFIELD .

Best regards, Kostya

[gevero]

Hi Kostya

Thanks for your comment. I must confess that a lot of time passed since I implemented the single particle routines of py_gmm, almost ten years, so I think I am kind of rusty... My understanding is that the problem is of course mathematical, but the stability (or instability) of the recurrence is determined by the fact that riccati-bessel functions of the first kind monotonically decrease as the order n increases. So, for upward recurrence, at each jump of order of magnitude, I am losing a significant digit. On the other hand, if I perform downward recurrence, even if starting from a random value, I gain a significant digit each jump of order of magnitude. Of course it may take many iterations to recover the needed significant digits. The paper I linked uses Kapteyn inequality to estimate the number of recursions needed to recover all the necessary significant digits, so It was my belief that py_gmm could handle the situation, nevertheless I might be wrong, especially because I almost always work with metals. In any case, I think that a size parameter of 10 or 20 should not be an issue for any normal Mie code. I shall setup a toy model very soon, in order to check my code.

Best

Giovanni

[kostyfisik]

Dear Giovanni,

Thank you for the link, one more (looks to be useful too) I found in your code, so both this papers are mentioned in an issue for scattnalay about large size parameters. https://github.com/ovidiopr/scattnlay/issues/3

I added a proof-of-concept multiprecision support (early alpha) to Scattnlay, with C++ templates and Boost::Multiprecision lib it was rather strighforward, the code is in "mp" branch at GitHub now . I hope I will have some spare time later on to provide a detailed research in this area and to share the results...

I have check Scattnlay with FEM for large (D=1.5\lambda) particles, and we have addittional self check due to multilayer property (the outgoing field from the inner layer should be zero, if it is not - for sure there is a stability problem). Still it looks to me like the multi precision is the only reliable way to do the check for any size parameter\material properties... Very fast FEM for large 3D particles becomes too time consuming and memory hungry...

Best regards, Kostya

[kostyfisik]

Hi all,

If anyone cares - it looks that multi-precision in Scattnlay works fine now, I was told by @ovidiopr that he checked several examples from Du`s papers, MP fit well while double precision had failed.

[gevero]

Hi Kostya

Thanks for the update. I am happy that Scattnlay work fine. Still MIECPP from Du's paper is a code working in double precision, so i do not understand why another code would need multiple precision to do the same computation. I still believe that the whole issue can be solved by a correct implementation in double precision.

Best

Giovanni

[kostyfisik]

For sure you are right, Du`s algorithm seems to be a superior one. However, it is always nice to verify the result using several independent approaches.

[gevero]

Sure!!! I browsed the code from scattnlay and I guess that the problem is two-fold:

1) the recursion for bessel functions of complex argument is performed directly, and this exposes the code to an overflow problem when the imaginary part of the argument is big. If instead one performs the recursion over the ratio of bessel functions, as in Du's paper, the problem is solved. One then of course needs to express the all expansion coefficients (a[n,l],b[n,l],...) for each layer of the sphere as a function of the ratio of bessel functions. 2) The recursion must be performed using kaptein inequality as in Du, in order to recover all the needed significant digits

Other than that no futher modifications should be needed

Best

Giovanni

[kostyfisik]

Giovanni! You are great! I hope I will find some free time to read Du paper carefully and to take you notes into account. However, this is not a plan for near future :(

[gevero]

I wish I was more familiar than I am with C++, otherwise I could try a pull request myself. Anyway never say never! :)

[kostyfisik]

Feel free to ask any help if you decide to start :) Actually, there are many C++ tutorials around, and modern C++11 with std:: library is very expressive, most of the functions needed for scientific computing you can find in std:: or boost:: libs, and there are a number of good other libs easily available under linux.

The most obvious benefit of C++ is performance - in some case you can turn your C++ to use all cores with an additional single line in the code (via OpenMP routines). With python it turns to be much harder.

[jeffmoreland]

Giovanni,

Did you ever create the toy example for Paul with large silica particles in water consisting of "very large radii (e.g. 5-20 wavelengths) and a refractive index of glass (~1.46) embedded in a medium (~1.333)"? I would be interested to see the code for that example.

[gevero]

Hi Jeff

Actually, I did not. But since you asked, I will post it so in the next couple of days.

Best

Giovanni