Please help with usage of fieldnlay
Closed this issue · 13 comments
I am having difficulties reproducing the phase retardation computed for a simple sphere. I would like to obtain the field behind a dielectric sphere with a radius of 5 wavelengths at a distance of 6 wavelengths from the center of the sphere.
Using the examples directory, I came up with an example script. Could you please help me to find out what is wrong with it? The retrieved phase retardation is too low (-0.28 to -0.21 rad) when compared to other methods like GMM-FIELD (0 to 0.64 rad).
import numpy as np
import scattnlay
import matplotlib.pylab as plt
n1 = 1.01 #weak dielectric object
nm = 1.0 #vacuum
radius = 5.0 #radius of the sphere in wavelengths?
extent = 20 #extent of the simulation size in wavelengths?
distance = 6.0 #distance where we want to have the measured field behind the sphere in wavelengths measured from the center of the sphere?
resolution = 10.0 #pixels per wavelength in the output image
# There is only one sphere, no layers
x = np.ones((1, 1), dtype = np.float64)
x[0, 0] = radius
# Set the refractive index of the sphere, normalized to that of the medium
m = np.ones((1, 1), dtype = np.complex128)
m[0, 0] = n1/nm
nptsx = extent*resolution
nptsy = extent*resolution
scanx = np.linspace(-extent/2, extent/2, nptsx, endpoint=True)
scany = np.linspace(-extent/2, extent/2, nptsy, endpoint=True)
coordX, coordY = np.meshgrid(scanx, scany)
coordX.resize(nptsx*nptsy)
coordY.resize(nptsx*nptsy)
coordZ = np.ones(nptsx*nptsy, dtype=np.float64)*distance
coord = np.vstack((coordX, coordY, coordZ)).transpose()
terms, E, H = scattnlay.fieldnlay(x, m, coord)
# take the x-component of the electric field
Ex = E[:,:,0].reshape(nptsx, nptsy)
# plot the phase of the x-component of the electric field
ax = plt.subplot(111)
mapper = plt.imshow(np.angle(Ex))
plt.colorbar(mapper, ax=ax, label="phase [rad]")
plt.title("scattering by a dielectric sphere")
plt.show()
# Output image attached
There seems to be a negative offset in the image. Maybe I need to perform some kind of background correction?
There is also a little confusion on my side about the variable x. In one example (https://github.com/ovidiopr/scattnlay/blob/master/examples/field-SiAgSi.py#L72-74), it us multiplied by a factor of 2PI. In another example (https://github.com/ovidiopr/scattnlay/blob/master/examples/field-silica.py#L43) it is not.
Thanks for your help,
Paul
Hi, Paul! Variable x is given in dimensionless units (size parameter). So you can set it directly or convert it from wavelength and particle size, given e.g. in nanometers. This is the difference in the examples comes from.
Talking about the main part of your question - in general, Mie codes (including GMM) do not behave quite well for large particles (I have verified Scattnlay with FEM and FDTD for size of the particle of up to 1.5 WL). So for 6 WL you can probably have some problems (due to the need of calculations of spherical harmonics), and you will have them for sure for realistic metals (with large imaginary part of index). In Scattnlay we used and advanced method for calculations for N and M via Ricatti-Bessel functions (Dn in the notation at https://github.com/ovidiopr/scattnlay/blob/master/src/nmie.cc#L728 ) and GMM uses Bessel and Hankel funcitons directly (see vswf.f in GMM source), so I would expect our code to give a more reliable result for large particles. Scattnlay has a built-in stability check at https://github.com/ovidiopr/scattnlay/blob/master/src/nmie.cc#L1132 (uncoment two throw lines below to be sure not to get unstable result), so check the output (I am not sure how does it work with python). The stability check simply comes from the fact that we have additional equations for expansion coefficients for outcoming wave in the inner layer (a_0 = 0 and b_0 = 0).
For bulk and core-shell cases you can also may want to check BHFIELD code by Suzuki http://www.scattport.org/index.php/light-scattering-software/mie-type-codes/list/521-bhfield It has an arbitraray precision mode (so you can change the number of digits used for calculations and check the stability this way). So Scattnlay and GMM uses double precision numbers (16 digits), with ARP version of BHFIELD you can check 50 and 100 digits precision to be sure with the convergence.
Best regards,
Kostya
#3 can be related
Thanks Kostya,
when I uncomment the two lines and run make python_ext again, nothing is thrown with the above example. Here is the result compared to that of GMM-FIELD. I still believe that there is something wrong with my script.
Regarding the 2PI: Is x supposed to be in wavelengths? If not, how is the wavelength defined? If yes, why is there a 2PI in the example, even though the radii are normalized to the wavelength WL.
Thank you again for clarifying. I might have a look at other codes. However, I would like to stay as close to Python as possible.
Kind Regards,
Paul
x - is a size parameter (is always used in all Mie books and papers). To convert from real world values you can use x=2_PI__R/WL where radius and wavelength are given in same units. R/WL = radius from your code, so you are probably have missed 2pi.
Hi Paul,
I have checked your code and Kostya is right, you are missing a 2PI. The relevant changes are:
radius = 2.0*np.pi*5.0 #radius of the sphere in wavelengths?
extent = 2.0*np.pi*20.0 #extent of the simulation size in wavelengths?
distance = 2.0*np.pi*6.0 #distance where we want to have the measured field behind the sphere in wavelengths measured from the center of the sphere?
Unfortunately the calculation becomes unstable for this combination of parameters (the program will give you a warning). Do you have an example with a lower size parameter (x) to verify your script?
Best regards,
Ovidio
Thank you both! That solved it. With the 2PI the program also aborts because of unstable calculations:
terminate called after throwing an instance of 'std::invalid_argument'
what(): Unstable calculation of aln_[0][n]!
Aborted (core dumped)
When I change to stable variables, scattnlay agrees very well with GMM-FIELD:
FYI: This is what happens when the radius of the sphere is increase to 5 wavelengths:
radius = 5
extent = 10
distance = 6
resolution = 5.0 
I am planning to look into bhfield to see where things go wrong.
This is the full example. I added a background correction step:
from __future__ import division
import numpy as np
import scattnlay
import matplotlib.pylab as plt
# weak dielectric object
n1 = 1.01
#surrounding medium is vacuum
nm = 1.0
# radius of the sphere in wavelengths
radius = .4
# extent of the simulation size in wavelengths
extent = 4
# distance where we want to have the measured field behind the sphere
# in wavelengths measured from the center of the sphere
distance = 0.4
# pixels per wavelength in the output image
resolution = 10.0
# size parameters need to be multiplied by 2 pi
twopi = 2*np.pi
# There is only one sphere, no layers
x = np.ones((1, 1), dtype = np.float64)
x[0, 0] = radius*twopi
# Set the refractive index of the sphere, normalized to that of the medium
m = np.ones((1, 1), dtype = np.complex128)
m[0, 0] = n1/nm
nptsx = extent*resolution
nptsy = extent*resolution
scanx = np.linspace(-extent/2, extent/2, nptsx, endpoint=True)*twopi
scany = np.linspace(-extent/2, extent/2, nptsy, endpoint=True)*twopi
coordX, coordY = np.meshgrid(scanx, scany)
coordX.resize(nptsx*nptsy)
coordY.resize(nptsx*nptsy)
coordZ = np.ones(nptsx*nptsy, dtype=np.float64)*distance*twopi
coord = np.vstack((coordX, coordY, coordZ)).transpose()
terms, E, H = scattnlay.fieldnlay(x, m, coord)
# take the x-component of the electric field
Ex = E[:,:,0].reshape(nptsx, nptsy)
# background
Ex /= np.exp(1j*2*np.pi*distance*nm)
# plot the phase of the x-component of the electric field
ax = plt.subplot(111)
mapper = plt.imshow(np.angle(Ex))
plt.colorbar(mapper, ax=ax, label="phase [rad]")
plt.title("scattering by a dielectric sphere")
plt.show()It would be cool, if you added that to the repository.
Cheers,
Paul
Paul, please, add copyright and license (we use GPLv3+) info and provide a pull request. It would be nice if you can avoid imports from future (still it runs well at my laptop so not a big problem). It would be ideal if you can provide few words about the problem. Actually, I do not understand the physical meaning of the plotted value.
Hello Kostya,
yes I will do all that. I will also choose values that are closer to an experimental situation.
I have another question regarding my example. When the refractive index of the medium is not 1, then nm does also go into the variable twopi above, right? i.e. twopi = 2*np.pi*nm
Here is a comparison with b/w gmm-field, scattnlay and bhfield:
Yes, media refractive index is also a part of a size parameter. You can check the attached paper for details, the evaluations in the Scattnlay are very close to the printed equations. How many digits you have used for BHFIELD? As far as it is exactly the same result as GMM it looks like you have messed versions and were using double precision BHFIELD instead of it`s ARP version.
Thanks.
I am running gmmfield at the maximum precision that (I think) can be set manually during compile time, i.e. in the file gmm01f.par I write parameter(nLp=2,np=157). I find that so far, gmm-field and bhfield agree perfectly in the small interval that I am testing in. Just for testing, I have set the parameter mpdigits to 200 for bhfield. The result is the same figure as above.
This is the version of gmmfield that I am using:
http://moritz-ringler.name/dissertation/GMM_FIELD.tar.bz2
More precisely, the refractive index can be part of the size parameter. It depends on the meaning of the wavelength for you. If you define lambda as the wavelength in vacuo then you must use the refractive index to calculate the size parameter. However, if you define lambda as the wavelength inside the medium then you don't have to use the refractive index to calculate the size parameter. The important thing is to be consistent once you have made a choice.
Regards,
Ovidio
Thanks. Yes, I tend to use the wavelength in vacuo as a unit length. I have updated the example to avoid any confusion in the future.
Kind Regards,
Paul
Paul, I have added multi precision support to Scattnlay, it still very early alpha, however, you are welcome to participate. I am using Boost::Multiprecision (so you need to install Boost libs to compile and run), the code is available at mp branch at GitHub. There are a lot of things to do in C++ part, I hope we can deal it a little bit later (like put Pi constant with needed precision, precise conversion of input and so on). I am always a bit confused with python part, so I was not able to make it run fast. Still is should not be hard as far as the API for C functions (used with Cython) is still the same.
To switch back to double of float precision simply change namespace back to std:: https://github.com/ovidiopr/scattnlay/blob/mp/src/nmie-impl.hpp#L59
and select float type
https://github.com/ovidiopr/scattnlay/blob/mp/src/nmie.cc#L54
(I have renamed header to *.hpp extension to keep C++ notation and to stress the fact that it is fully template based now).
Test program output for 100 digits of precision...
~/symlink/scattnlay/tests/c++$ ./go-speed-test.sh
Compile with gcc
Should be:
test01, +1.41154e+00, +4.17695e-01, +9.93844e-01, +1.59427e-01, +1.25809e+00, +3.67376e-01, +2.95915e-01
Result:
Qext = 1.411538881495128733962719194753380836603011775206276379868166282271624619065332987629814088255337614
-- C++ time consumed 0.0147995 sec
test01, +1.41154e+00, +4.17695e-01, +9.93844e-01, +1.59427e-01, +1.25809e+00, +3.67376e-01, +2.95915e-01
Theta, S1.r, S1.i, S2.r, S2.i
0.00, +3.52885e-01, -4.27839e-01, +3.52885e-01, -4.27839e-01
22.50, +3.44554e-01, -4.14362e-01, +3.29828e-01, -3.94049e-01
45.00, +3.21195e-01, -3.76903e-01, +2.65943e-01, -3.03992e-01
67.50, +2.87219e-01, -3.23325e-01, +1.75076e-01, -1.85579e-01
90.00, +2.48609e-01, -2.63824e-01, +7.48392e-02, -6.95864e-02