ovidiopr/scattnlay

Near- and far-field Mie scattering by a multilayered sphere

scattnlay's Issues

Large complex argument

At the moment Scattnlay is limited in the range of used materials due to constraints on accuracy of calculation of Bessel functions with large complex argument. May be we can rewrite code to benefit from using psi_(n-1)/psi_n instead of diff(psi)/psi to calculate scattering coefficients with large complex argument? See H. Du, "Mie-scattering calculation," Appl. Opt. 43, 1951-1956 (2004). http://dx.doi.org/10.1364/AO.43.001951 http://www.opticsinfobase.org/ao/abstract.cfm?uri=ao-43-9-1951

Please help with usage of fieldnlay

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

errors when trying to install scattnlay in Mac

In Mac M2 Ventura 13.6.1 (Xcode version is 15.2), I got many errors when executing in the terminal:

python3 -m pip install python-scattnlay-2.3.tar.gz

/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/c++/v1/cmath:607:12: error: no matching function for call to 'isinf'

/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/c++/v1/cmath:587:12: error: no matching function for call to 'isnan'

/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/c++/v1/complex:637:19: error: no matching function for call to '__constexpr_copysign'

/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/c++/v1/complex:638:19: error: no matching function for call to '__constexpr_copysign'

/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/c++/v1/complex:706:63: error: no matching function for call to '__constexpr_fabs'

/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/c++/v1/cmath:627:31: error: floating point classification requires argument of floating point type (passed in 'boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<100>>')

/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/c++/v1/complex:762:19: error: no matching function for call to '__constexpr_copysign'

fatal error: too many errors emitted, stopping now [-ferror-limit=]
2 warnings and 20 errors generated.
error: command '/usr/bin/clang' failed with exit code 1
[end of output]

note: This error originates from a subprocess, and is likely not a problem with pip.
ERROR: Failed building wheel for python-scattnlay

switching from Cython to pybind11?

Additional tests

There is also a number of stupid tests which should be easily passed with the correct code, e.g.: several layers with the same index should give the same result for the same outer size, adding one more outer layer with air index should not change the result, tangential field components to the left and to the right of the border between layer should be the same, if we decrease some layer thickness to zero (or to almost zero) the changes of the index value of this layer should not change the result and so on...

Microsphere nanojet simulation

Hi! How can I increase the nmax in python code if I compute a large-scale (~ 10 wavelengths )sphere ? I am looking forward to your reply! Thank you!

Improvements for webapp

• Far-field plot + polarization options
• (bug) For multiple spectrum plots wrong conversion of plane wave unit.
• (bug) GaP_Khmelevskaia - wrong shown range in material explorer?
• Near-field - split save button for data and plotting script
• Near-field (IR) - fixed -> fixed Lx/Ly
• Near-field (IR) - put EHk vectors to the plot
• Near-field (IR) - Find a better name for 'Plot relative center'
• Near-field (IR) - separate "particle setup" from "spectrum" to a separate tab.
• Near-field (IR) - show number of multipoles used for evaluation
• Near-field - power-flow lines - normalize vector (we just need the direction) and use boost::numeric::odeint using direction vector as a derivative (so it will interpret vector direction as a unit speed vector).
• Near-field - show field without incident wave
• Add search with logical AND to the q-table with all materials
• Vuepress with Vite documentation for webapp
• Switch between basic and advanced modes
• Near-field - 3D vector plot, e.g. https://en.wikipedia.org/wiki/Vector_spherical_harmonics#/media/File:VSHwiki.svg
• Near-field - vector streamline plot for the field
• Evaluate poles of scattering matrix (KK) - can be a problem with dispersion
• Near-field - multipole decomposition
• Near-field - magnetic field
• Use BEM for CSS styles (CSS modules?)
• Set nmax for simulation.
• Set dielectric material with eps and tan(delta)
• Verify user input (number of mode should be a positive integer, fromWL < toWL, etc.)
• Notify "Scroll down to see the result" after spectrum simulation
• Use https://troisjs.github.io/ for far-field visualization? Or D3FC WebGL for heatmap (near-field) for fast zoom and colorscale changes, e.g. https://blog.scottlogic.com/2020/05/01/rendering-one-million-points-with-d3.html
• Add presets to reproduce figures from papers
• Add cookies to remember user setting from previous visit + reset button.
• Add promotional video with nice girls using Mie calculator
• Solve inverse problem (see PyMieScatt for examples)
• Add magnetic materials
• Near-field - add save data button
• Near-field - additional settings to specify simulation plot area
• Switch between frequency/wavelength
• Add only n data, and fused silicon dioxide to the https://refractiveindex.info/?shelf=main&book=SiO2&page=Malitson
• Add button to view (as a table) and download (as a file) simulation results data.
• Add dispersion
• Near-field plot
• Add host media
• Additional WLs array for plotting (so changing source units should not be effect change of from\to\step WL from previous simulation)

Near-field verification for a 10mkm water droplet - create test case

Create test data with mpmath and set GTest for the test case described below (see #42 for more context)

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.

build errors, nscattlay_dp not found

System information:

Ubuntu 19.10
GCC 8.3.0 (installed with apt)
Boost 1.6.7 (installed with apt)
scattnlay - from github repo
Python 3.7 using miniconda

make ext  # works
make install # fails 

error message from make install:

`which python3` setup.py install --root / 
running install
running build
running build_py
error: package directory 'scattnlay_dp' does not exist
make: *** [Makefile:27: install] Error 1

Here is the directory output after make ext:

(base) vince@vince-vb:~/Coding/scattnlay$ ls -la
total 15632
drwxr-xr-x 11 vince vince    4096 Jun 24 13:10 .
drwxr-xr-x  6 vince vince    4096 Jun 24 11:24 ..
drwxr-xr-x  4 vince vince    4096 Jun 24 13:06 build
-rw-r--r--  1 vince vince       0 Jun 24 11:24 CHANGES
-rw-r--r--  1 vince vince    1404 Jun 24 11:24 COPYING
drwxr-xr-x  3 vince vince    4096 Jun 24 11:24 debian
drwxr-xr-x  2 vince vince    4096 Jun 24 11:24 doc
drwxr-xr-x  2 vince vince    4096 Jun 24 11:24 examples
drwxr-xr-x  8 vince vince    4096 Jun 24 11:24 .git
-rw-r--r--  1 vince vince     483 Jun 24 11:24 .gitignore
-rwxr-xr-x  1 vince vince    4528 Jun 24 11:24 go.sh
-rw-r--r--  1 vince vince     922 Jun 24 13:10 install.md
-rw-r--r--  1 vince vince   35147 Jun 24 11:24 LICENSE
-rw-r--r--  1 vince vince    2522 Jun 24 11:24 Makefile
-rw-r--r--  1 vince vince     284 Jun 24 11:24 MANIFEST.in
-rw-r--r--  1 vince vince     864 Jun 24 11:24 PKG-INFO
-rwxr-xr-x  1 vince vince     137 Jun 24 11:24 push-to-github.sh
-rw-r--r--  1 vince vince    5992 Jun 24 11:24 README.md
drwxr-xr-x  2 vince vince    4096 Jun 24 12:47 scattnlay
-rwxr-xr-x  1 vince vince 6494096 Jun 24 13:06 scattnlay_dp.cpython-37m-x86_64-linux-gnu.so
-rwxr-xr-x  1 vince vince 9373200 Jun 24 13:07 scattnlay_mp.cpython-37m-x86_64-linux-gnu.so
-rw-r--r--  1 vince vince      40 Jun 24 11:24 setup.cfg
-rwxr-xr-x  1 vince vince    3711 Jun 24 11:24 setup.py
drwxr-xr-x  2 vince vince    4096 Jun 24 11:24 src
drwxr-xr-x  5 vince vince    4096 Jun 24 11:24 tests
drwxr-xr-x  3 vince vince    4096 Jun 24 11:24 utils

Add support of azimutal and radial polarization of incident light.

E.g. try to reproduce results from https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.241110

Jump of field value at z=0

Consider the following command:

$ ./fieldnlay-mp -l 1 113 1.3346 0 -p 242.1 242.1 1 -187.1 -187.1 1 -0.001 0.001 9

EM field in the range of z∈[-0.001,0.001] should be very smooth, but the calculated values exhibit a jump at z=0:

Econv:5.62569 Hconv:1.77823
Econv:0.00969625 Hconv:0.0109756
Econv:0.00457643 Hconv:0.00520441
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
242.1000000, -187.1000000, -0.0010000, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan
242.1000000, -187.1000000, -0.0007500, +2.40842e-02, -1.73683e-02, -7.43300e-03, -2.51150e-02, -6.98682e-03, -2.42827e-02, +2.85789e-05, +2.27491e-05, +9.94119e-05, -1.28891e-05, -3.36413e-05, -3.00095e-05
242.1000000, -187.1000000, -0.0005000, +2.40856e-02, -1.73637e-02, -7.43103e-03, -2.51140e-02, -6.98698e-03, -2.42823e-02, +2.85769e-05, +2.27549e-05, +9.94121e-05, -1.28738e-05, -3.36341e-05, -3.00127e-05
242.1000000, -187.1000000, -0.0002500, +2.40870e-02, -1.73591e-02, -7.42906e-03, -2.51130e-02, -6.98714e-03, -2.42819e-02, +2.85749e-05, +2.27608e-05, +9.94124e-05, -1.28585e-05, -3.36269e-05, -3.00160e-05
242.1000000, -187.1000000,  0.0000000, +5.41821e-03, -1.86391e-02, +7.00168e-03, -2.41193e-02, -6.98730e-03, -2.42815e-02, +5.87756e-05, +2.29790e-06, +7.60713e-05, +2.97548e-06, -3.36197e-05, -3.00193e-05
242.1000000, -187.1000000,  0.0002500, +2.40899e-02, -1.73500e-02, -7.42513e-03, -2.51110e-02, -6.98746e-03, -2.42811e-02, +2.85708e-05, +2.27725e-05, +9.94129e-05, -1.28279e-05, -3.36124e-05, -3.00226e-05
242.1000000, -187.1000000,  0.0005000, +2.40913e-02, -1.73454e-02, -7.42316e-03, -2.51100e-02, -6.98762e-03, -2.42808e-02, +2.85688e-05, +2.27783e-05, +9.94131e-05, -1.28126e-05, -3.36052e-05, -3.00258e-05
242.1000000, -187.1000000,  0.0007500, +2.40928e-02, -1.73408e-02, -7.42120e-03, -2.51090e-02, -6.98777e-03, -2.42804e-02, +2.85667e-05, +2.27842e-05, +9.94134e-05, -1.27973e-05, -3.35980e-05, -3.00291e-05
242.1000000, -187.1000000,  0.0010000, +2.40942e-02, -1.73363e-02, -7.41923e-03, -2.51080e-02, -6.98793e-03, -2.42800e-02, +2.85647e-05, +2.27900e-05, +9.94136e-05, -1.27819e-05, -3.35908e-05, -3.00324e-05

Here's a plot of these values (51 points instead of the 9 above) to make the jump more obvious:

Support of materials with negative refractive index

To do:

• Verify math
• Add and verify example

size parameter should be a complex variable to deal with lossy host media.

Inconsistent results between PyPi package, python package installed from repository and web app

Dear devs,
i was having some trouble reproducing some Mie results with you PyPi package. After some digging i found an inconsistency between the PyPi package, the python package installed from the github repository and the web app.

minimal examples:
a) a sphere with static n=1.4727 + 0.00001j with radius 400nm at wavelength 800nm:

• same result using your web app and the python package install from the repository (version 2.3)
• different result using the PyPi installation (also version 2.3)

b) same sphere with static (n=1.4727 + 0.00001j with radius 400nm) at wavelength 8000 nm:

• all 3 methods give the same result

I used the exact same script with the python packages (PyPi vs GitHub) so I assume there is a bug somewhere in the PyPi package .

Christian

Error in the radiative pressure efficiency factor

Hello everybody,

I have seen something seeming like an error, in scattnlay v1.0.0. It has stayed in the newest codes.
The radiative pressure efficiency factor Qpr is finally calculated using the Eq. 29 of the Pena & Pal paper. At a point in the formula appears the ratio n*(n+2)/(n+1) .
In the code this is written with ... (n*(n + 2)/(n + 1)) ... , where n is an integer variable. The C/C++ compiler places here integer operations, including an integer division.
Such an operation discards the fractional part, and reading the publication it does not seem to be the wanted behaviour (since it is not mentioned in the paper formula, I suppose the division is an ordinary one). It would not be a problem if n*(n+2)/(n+1) had a null fractional part, but as n is a strictly positive integer this fractional part is never zero.
I tested that point a little bit, experiencing a difference of a few percents in Qpr when enforcing or not a floating point division.

Best regards,

memory required

Is it possible to estimate the amount of memory required when using standalone program fieldnlay.exe?

Update convergence preset?

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–7229 (2014).

nan values in the central part of the sphere

Hello.
I am using Windows 10. I have compiled scattnlay and I use standalone program fieldnlay.exe to calculate the electric field inside the sphere.
I use the command:
fieldnlay.exe -l 1 250 1.1 0 -p 0 0 1 -250 250 800 -250 250 800.
Calculations are successful but I have several elements with nan values in center of the sphere. The output table is as follows (for brevity, I do not list all the lines with nan elements):
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
...
-17.5000000, -2.5000000, -2.5000000, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan
-17.5000000, -2.5000000, 2.5000000, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan
-17.5000000, 2.5000000, -2.5000000, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan, +nan
...
Could you help me with this issue.
Thank you in advance.

Command line input

Shoul we use something standalone version for parameters intput? Like http://sourceforge.net/projects/tclap/

possible numerical problems

I am trying to use scattnlay, to reproduce figure 4.6 from Bohren & Huffman looking at the scattering properties of water droplets.

I pulled optical properties of water from Wikipedia:
https://en.wikipedia.org/wiki/Optical_properties_of_water_and_ice

There appear to be some numerical errors in calculating the Qext. I am wondering if this is common / expected or if it is a bug or if I am doing something wrong:

The following code was run in a Jupyter Notebook:

%pylab inline

import sys
sys.path.append('/home/vince/Coding/scattnlay/')
import scattnlay

wl_um = \
array([0.2  , 0.225, 0.25 , 0.275, 0.3  , 0.325, 0.35 , 0.375, 0.4  ,
       0.425, 0.45 , 0.475, 0.5  , 0.525, 0.55 , 0.575, 0.6  , 0.625,
       0.65 , 0.675, 0.7  , 0.725, 0.75 , 0.775, 0.8  , 0.825, 0.85 ,
       0.875, 0.9  , 0.925, 0.95 , 0.975, 1.   , 1.2  , 1.4  , 1.6  ,
       1.8  , 2.   , 2.2  , 2.4  , 2.6  , 2.65 , 2.7  , 2.75 , 2.8  ,
       2.85 , 2.9  , 2.95 , 3.   , 3.05 , 3.1  , 3.15 , 3.2  , 3.25 ,
       3.3  , 3.35 , 3.4  , 3.45 , 3.5  , 3.6  , 3.7  , 3.8  , 3.9  ,
       4.   , 4.1  , 4.2  , 4.3  , 4.4  , 4.5  , 4.6  , 4.7  , 4.8  ,
       4.9  , 5.   , 5.1  , 5.2  , 5.3  , 5.4  , 5.5  , 5.6  , 5.7  ,
       5.8  , 5.9  , 6.   , 6.1  , 6.2  , 6.3  , 6.4  , 6.5  , 6.6  ,
       6.7  , 6.8  , 6.9  , 7.   ])

n_water = \
array([1.396, 1.373, 1.362, 1.354, 1.349, 1.346, 1.343, 1.341, 1.339,
       1.338, 1.337, 1.336, 1.335, 1.334, 1.333, 1.333, 1.332, 1.332,
       1.331, 1.331, 1.331, 1.33 , 1.33 , 1.33 , 1.329, 1.329, 1.329,
       1.328, 1.328, 1.328, 1.327, 1.327, 1.327, 1.324, 1.321, 1.317,
       1.312, 1.306, 1.296, 1.279, 1.242, 1.219, 1.188, 1.157, 1.142,
       1.149, 1.201, 1.292, 1.371, 1.426, 1.467, 1.483, 1.478, 1.467,
       1.45 , 1.432, 1.42 , 1.41 , 1.4  , 1.385, 1.374, 1.364, 1.357,
       1.351, 1.346, 1.342, 1.338, 1.334, 1.332, 1.33 , 1.33 , 1.33 ,
       1.328, 1.325, 1.322, 1.317, 1.312, 1.305, 1.298, 1.289, 1.277,
       1.262, 1.248, 1.265, 1.319, 1.363, 1.357, 1.347, 1.339, 1.334,
       1.329, 1.324, 1.321, 1.317])

k_water = \
array([1.10e-07, 4.90e-08, 3.35e-08, 2.35e-08, 1.60e-08, 1.08e-08,
       6.50e-09, 3.50e-09, 1.86e-09, 1.30e-09, 1.02e-09, 9.35e-10,
       1.00e-09, 1.32e-09, 1.96e-09, 3.60e-09, 1.09e-08, 1.39e-08,
       1.64e-08, 2.23e-08, 3.35e-08, 9.15e-08, 1.56e-07, 1.48e-07,
       1.25e-07, 1.82e-07, 2.93e-07, 3.91e-07, 4.86e-07, 1.06e-06,
       2.93e-06, 3.48e-06, 2.89e-06, 9.89e-06, 1.38e-04, 8.55e-05,
       1.15e-04, 1.10e-03, 2.89e-04, 9.56e-04, 3.17e-03, 6.70e-05,
       1.90e-02, 5.90e-02, 1.15e-01, 1.85e-01, 2.68e-01, 2.98e-01,
       2.72e-01, 2.40e-01, 1.92e-01, 1.35e-01, 9.24e-02, 6.10e-02,
       3.68e-02, 2.61e-02, 1.95e-02, 1.32e-02, 9.40e-03, 5.15e-03,
       3.60e-03, 3.40e-03, 3.80e-03, 4.60e-03, 5.62e-03, 6.88e-03,
       8.45e-03, 1.03e-02, 1.34e-02, 1.47e-02, 1.57e-02, 1.50e-02,
       1.37e-02, 1.24e-02, 1.11e-02, 1.01e-02, 9.80e-03, 1.03e-02,
       1.16e-02, 1.42e-02, 2.03e-02, 3.30e-02, 6.22e-02, 1.07e-01,
       1.31e-01, 8.80e-02, 5.70e-02, 4.49e-02, 3.92e-02, 3.56e-02,
       3.37e-02, 3.27e-02, 3.22e-02, 3.20e-02])

xvals = 1 / wl_um
fig, (ax_n, ax_k) = subplots(2, 1, figsize=(6,6))
ax_n.plot(xvals, n_water)
ax_k.plot(xvals, k_water, 'r')
fig.tight_layout()
ax_k.set_xlabel('inverse microns')
ax_k.set_ylabel('k')
ax_n.set_ylabel('n')
ax_n.set_title('refractive index of water')

radius = 1.0  # um
N_medium = 1  # air

x_water = (2 * pi * N_medium * radius / wl_um).reshape(-1, 1)
m = ((n_water + 1j * k_water) / N_medium).reshape(-1, 1)
results = scattnlay.scattnlay(x_water, m)
Qext = results[1]

# create a finer grid for interpolation of refractive index properties
wl_fine = linspace(wl_um.min(), wl_um.max(), 1000)
n_water_fine = interp(wl_fine, wl_um, n_water)
k_water_fine = interp(wl_fine, wl_um, k_water)

x_water_fine = (2 * pi * N_medium * radius / wl_fine).reshape(-1, 1)
m_fine = ((n_water_fine + 1j * k_water_fine) / N_medium).reshape(-1, 1)
results_fine = scattnlay.scattnlay(x_water_fine, m_fine)
Qext_fine = results_fine[1]

fig, (ax1, ax2) = subplots(1, 2, figsize=(9, 4), dpi=90, facecolor='white')
ax1.plot(1/wl_fine, Qext_fine, label='interpolated')
ax1.plot(1/wl_um, Qext, label='table')
ax2.semilogy(1/wl_fine, Qext_fine, label='interpolated')
ax2.semilogy(1/wl_um, Qext, 'o', label='table')
for ax in (ax1, ax2):
    ax.set_xlabel('inverse microns')
    ax.set_ylabel('Excinction efficiency')
    ax.legend()
fig.tight_layout()

idx = argmin(Qext)
ax1.text(1/wl_um[idx], Qext[idx],
         f"neg Qext @ {wl_um[idx]} microns",
         horizontalalignment='left')

Fused silicon dioxide

Hi!
Could you add Fused silicon dioxide to the default list? the url from the refractiveindex.info is https://refractiveindex.info/?shelf=main&book=SiO2&page=Malitson

the simulator doesn´t let me add it
thank you!!

First point calculated by fieldnlay is always inf&nan

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.

Compilation issues for mac

Hi,

I am trying to compile the code for python in Mac OS 14.2.1. I have installed boost c++ library as well and directed to the include directory for the resource files. At least from github it also says that all checks have failed as well.

• When running pip install --global-option=build_ext --global-option="-I/opt/homebrew/Cellar/boost/1.83.0/include/" . --user I encounter this error:

/Users/joebiv98/anaconda3/lib/python3.10/site-packages/setuptools/distutils/dist.py:265: UserWarning: Unknown distribution option: 'extra_compile_args'
warnings.warn(msg)
In file included from src/nmie.cc:40:
src/nmie.h:150:27: warning: default member initializer for non-static data member is a C++11 extension [-Wc++11-extensions]
bool isExpCoeffsCalc = false;
^
src/nmie.h:151:27: warning: default member initializer for non-static data member is a C++11 extension [-Wc++11-extensions]
bool isScaCoeffsCalc_ = false;
^
src/nmie.h:152:27: warning: default member initializer for non-static data member is a C++11 extension [-Wc++11-extensions]
bool isMieCalculated_ = false;
^
src/nmie.h:156:29: warning: default member initializer for non-static data member is a C++11 extension [-Wc++11-extensions]
int PEC_layer_position_ = -1;
^
src/nmie.h:159:15: warning: default member initializer for non-static data member is a C++11 extension [-Wc++11-extensions]
int nmax_ = -1;
^
src/nmie.h:160:22: warning: default member initializer for non-static data member is a C++11 extension [-Wc++11-extensions]
int nmax_preset_ = -1;
^
src/nmie.h:166:18: warning: default member initializer for non-static data member is a C++11 extension [-Wc++11-extensions]
double Qsca_ = 0.0, Qext_ = 0.0, Qabs_ = 0.0, Qbk_ = 0.0, Qpr_ = 0.0, asymmetry_factor_ = 0.0, albedo_ = 0.0;
^
src/nmie.h:166:31: warning: default member initializer for non-static data member is a C++11 extension [-Wc++11-extensions]
double Qsca_ = 0.0, Qext_ = 0.0, Qabs_ = 0.0, Qbk_ = 0.0, Qpr_ = 0.0, asymmetry_factor_ = 0.0, albedo_ = 0.0;
^
src/nmie.h:166:44: warning: default member initializer for non-static data member is a C++11 extension [-Wc++11-extensions]
double Qsca_ = 0.0, Qext_ = 0.0, Qabs_ = 0.0, Qbk_ = 0.0, Qpr_ = 0.0, asymmetry_factor_ = 0.0, albedo_ = 0.0;
^
src/nmie.h:166:56: warning: default member initializer for non-static data member is a C++11 extension [-Wc++11-extensions]
double Qsca_ = 0.0, Qext_ = 0.0, Qabs_ = 0.0, Qbk_ = 0.0, Qpr_ = 0.0, asymmetry_factor_ = 0.0, albedo_ = 0.0;
^
src/nmie.h:166:68: warning: default member initializer for non-static data member is a C++11 extension [-Wc++11-extensions]
double Qsca_ = 0.0, Qext_ = 0.0, Qabs_ = 0.0, Qbk_ = 0.0, Qpr_ = 0.0, asymmetry_factor_ = 0.0, albedo_ = 0.0;
^
src/nmie.h:166:93: warning: default member initializer for non-static data member is a C++11 extension [-Wc++11-extensions]
double Qsca_ = 0.0, Qext_ = 0.0, Qabs_ = 0.0, Qbk_ = 0.0, Qpr_ = 0.0, asymmetry_factor_ = 0.0, albedo_ = 0.0;
^
src/nmie.h:166:108: warning: default member initializer for non-static data member is a C++11 extension [-Wc++11-extensions]
double Qsca_ = 0.0, Qext_ = 0.0, Qabs_ = 0.0, Qbk_ = 0.0, Qpr_ = 0.0, asymmetry_factor_ = 0.0, albedo_ = 0.0;
^
src/nmie.cc:293:10: warning: 'auto' type specifier is a C++11 extension [-Wc++11-extensions]
for (auto f:E)
^
src/nmie.cc:293:16: warning: range-based for loop is a C++11 extension [-Wc++11-extensions]
for (auto f:E)
^
src/nmie.cc:296:10: warning: 'auto' type specifier is a C++11 extension [-Wc++11-extensions]
for (auto f:H)
^
src/nmie.cc:296:16: warning: range-based for loop is a C++11 extension [-Wc++11-extensions]
for (auto f:H)
^
src/nmie.cc:304:29: error: expected expression
ml_mie.SetFieldCoords({Xp_vec, Yp_vec, Zp_vec});
^
src/nmie.cc:426:10: warning: 'auto' type specifier is a C++11 extension [-Wc++11-extensions]
for (auto curr_layer_size : layer_size) {
^
src/nmie.cc:426:31: warning: range-based for loop is a C++11 extension [-Wc++11-extensions]
for (auto curr_layer_size : layer_size) {
^
src/nmie.cc:1003:25: warning: comparison of integers of different signs: 'int' and 'size_type' (aka 'unsigned long') [-Wsign-compare]
for (int t = 0; t < theta_.size(); t++) {
~ ^ ~~~~~~~~~~~~~
src/nmie.cc:1080:5: warning: 'auto' type specifier is a C++11 extension [-Wc++11-extensions]
auto& m = refractive_index_;
^
src/nmie.cc:1168:40: error: non-aggregate type 'std::vector<std::complex >' cannot be initialized with an initializer list
std::vector<std::complex > ipow = {c_one, c_i, -c_one, -c_i}; // Vector containing precomputed integer powers of i to avoid computation
^ ~~~~~~~~~~~~~~~~~~~~~~~~~~
src/nmie.cc:1265:10: warning: 'auto' type specifier is a C++11 extension [-Wc++11-extensions]
for (auto& f : E_) f.resize(3);
^
src/nmie.cc:1265:18: warning: range-based for loop is a C++11 extension [-Wc++11-extensions]
for (auto& f : E_) f.resize(3);
^
src/nmie.cc:1266:10: warning: 'auto' type specifier is a C++11 extension [-Wc++11-extensions]
for (auto& f : H_) f.resize(3);
^
src/nmie.cc:1266:18: warning: range-based for loop is a C++11 extension [-Wc++11-extensions]
for (auto& f : H_) f.resize(3);
^
25 warnings and 2 errors generated.
error: command '/usr/bin/clang' failed with exit code 1
[end of output]

note: This error originates from a subprocess, and is likely not a problem with pip.
ERROR: Failed cleaning build dir for python-scattnlay
Failed to build python-scattnlay
ERROR: Could not build wheels for python-scattnlay, which is required to install pyproject.toml-based projects

My python version is 3.10.11, and my pip installer is 23.3.2. Any advice would be much appreciated! Thank you.

Compilation problem

Hello.
I want to compile code for using in python project in Windows 10. I work in msys command shell.

1. I downloaded scattnlay project.
2. Run python -m pip install numpy command. Output is 'Successfully installed numpy-1.17.2'.
3. Run make ext command.
   Output:
   '/c/Users/konstantin/AppData/Local/Programs/Python/Python37-32/python.exe' setup.py build_ext --inplace
   Traceback (most recent call last):
   File "setup.py", line 44, in
   import pybind11 as pb
   ModuleNotFoundError: No module named 'pybind11'
   make: *** [ext] Error 1
4. Run pip install pybind11 command. Output is 'Successfully installed pybind11-2.4.2'.
5. Run make ext command.
   Output:
   '/c/Users/konstantin/AppData/Local/Programs/Python/Python37-32/python.exe' setup.py build_ext --inplace
   running build_ext
   building 'scattnlay_dp' extension
   error: Unable to find vcvarsall.bat
   make: *** [ext] Error 1

Could you help me. Maybe I made mistakes in the earlier steps. And I have no idea about installing python-all-dev and python-numpy-dev. Is it enough to call a command python -m pip install numpy?

I also tried to install packege python-scattnlay (v. 2.2) in Pycharm project. But error occured (error: Microsoft Visual C++ 14.0 is required). Is it possible to avoid installation Microsoft Visual C++ 14.0 if I have gcc compiler.

Thank you in advance.

Add v2.0 tag for new relase

$ git tag -a v2.0 -m 'Scattnlay 2.0'

Try to use arblib

arblib can at least be used for creating test data
https://arblib.org/acb_hypgeom.html?highlight=bessel#c.acb_hypgeom_bessel_i

another option - user arblib to get a starting value of D1 for backward recurrence.

pip version tests fail with buffer wrong dimensions

System Information:

Ubuntu 19.10
GCC 8.3.0 (installed with apt)
Boost 1.6.7 (installed with apt)
scattnlay - installed via pip install
Python 3.7 using miniconda

Hi I am exploring this package, but am unable to run the tests. For example, when I run:

~/scattnlay/tests/python/>>> python test01.py

it prints the error message:

Traceback (most recent call last):
  File "test01.py", line 72, in <module>
    terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
  File "src/scattnlay.pyx", line 97, in scattnlay.scattnlay (src/scattnlay.cpp:2577)
ValueError: Buffer has wrong number of dimensions (expected 2, got 1)

I get a similar error when trying to run examples/calcspectra.py

Traceback (most recent call last):
  File "calc-spectra.py", line 62, in <module>
    mp=True
  File "../scattnlay/main.py", line 137, in scattnlay
    return scattnlay_(x, m, theta, nmax=nmax, pl=pl)
  File "src/scattnlay_mp.pyx", line 97, in scattnlay_mp.scattnlay (src/scattnlay_mp.cpp:2574)
ValueError: Buffer has wrong number of dimensions (expected 2, got 1)

DOI for scattnay

There is a possibility to apply for a DOI to cite GitHub release of Scattnlay in any futher publications. https://guides.github.com/activities/citable-code/

SetMaxTerms(-1) yields wrong results

I'm not sure what motivation for #23 was, but the following gives a test case that shows that current automatic determination of the number of required terms works incorrectly.

Consider the following calculation:

$ ./fieldnlay-mp -l 1 2 1.3 0 -p 1e-3 1e-3 1 1e-3 1e-3 1 -200 200 2000 | cut -d, -f3,4 > /tmp/field.csv

This gives us ℜ{Eₓ} vs z table. The result looks as follows:

This doesn't make much sense: we should see incident wave with a bit of interference due to the scattered wave, rather than the widely-changing envelope we actually got. Let's now explicitly set 100 terms in nearfield.cc:

diff --git a/src/nearfield.cc b/src/nearfield.cc
index e842115..a914248 100644
--- a/src/nearfield.cc
+++ b/src/nearfield.cc
@@ -235,7 +235,7 @@ int main(int argc, char *argv[]) {
       }
     }
 
-    nmie::nField(L, -1, x, m, -1, nmie::Modes::kAll, nmie::Modes::kAll, total_points, Xp, Yp, Zp, E, H);
+    nmie::nField(L, -1, x, m, 100, nmie::Modes::kAll, nmie::Modes::kAll, total_points, Xp, Yp, Zp, E, H);
 
     if (has_comment)
       printf("%6s\n", comment.c_str());

Now the results look much more sensible. The incident wave becomes visible much farther from the sphere, although it still goes wild at some distance:

Changing the maximum to 250 gives the good expected result: