merge

.gitignore

@@ -27,9 +27,18 @@
 *.out
 *.app
 
+fieldnlay-dp
+fieldnlay-mp
+scattnlay-dp
+scattnlay-mp
+
 # Binary
 *.bin
 
+# CLion files
+.idea
+cmake-build-*
+
 # Temp
 *~
 oprofiletmp
@@ -72,7 +81,6 @@ yarn-debug.log*
 yarn-error.log*
 
 # Editor directories and files
-.idea
 .vscode
 *.suo
 *.ntvs*
@@ -82,4 +90,3 @@ yarn-error.log*
 
 # end of web vue cli section
 ###########################################
-

CMakeLists.txt

@@ -0,0 +1,113 @@
+cmake_minimum_required(VERSION 3.15)
+project(scattnlay VERSION 2.3)
+
+cmake_host_system_information(RESULT HOSTNAME QUERY HOSTNAME)
+
+message("Build type is: ${CMAKE_BUILD_TYPE}")
+message("Host OS System: ${CMAKE_HOST_SYSTEM}")
+message("Hostname:  ${HOSTNAME}")
+message("CMake version:  ${CMAKE_VERSION}")
+
+# Select flags.
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O3 -g ")
+set(CMAKE_CXX_FLAGS_RELEASE "-O3")
+set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g")
+
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
+add_compile_options(-Wall)
+add_compile_options(-funroll-loops -fstrict-aliasing)
+
+# Set options
+option(ENABLE_MP "Use multiple precision" OFF)
+if (${ENABLE_MP})
+    add_compile_options(-DMULTI_PRECISION=100)
+endif ()
+
+# compiler details
+message("  C++ Compiler: ${CMAKE_CXX_COMPILER_ID} "
+        "${CMAKE_CXX_COMPILER_VERSION} "
+        "${CMAKE_CXX_COMPILER_WRAPPER}")
+
+# installation details
+message("  Installation prefix: ${CMAKE_INSTALL_PREFIX}")
+
+# Find Boost
+set(BOOSTROOT $ENV{BOOST_DIR})
+if (USE_STATIC_LIBRARIES)
+    set(Boost_USE_STATIC_LIBS ON)
+endif ()
+set(Boost_USE_MULTITHREADED OFF)
+set(Boost_USE_STATIC_RUNTIME OFF)
+find_package(Boost REQUIRED)
+if (Boost_INCLUDE_DIRS)
+    if (${Boost_VERSION} VERSION_LESS 1.60.0)
+        message(FATAL_ERROR
+                "Found Boost library is too old; required is version 1.60.0 or newer!")
+    endif ()
+    message("Found Boost include dir: ${Boost_INCLUDE_DIR}")
+    message("Found Boost library dir: ${Boost_LIBRARY_DIR}")
+    message("Found Boost libraries: ${Boost_LIBRARIES}")
+    include_directories(${Boost_INCLUDE_DIRS})
+endif ()
+
+#Find Python, NumPy and PyBind11
+find_package(Python COMPONENTS Interpreter Development)
+
+include_directories(${Python_INCLUDE_DIRS})
+
+message("Python_EXECUTABLE: ${Python_EXECUTABLE}")
+message("Python_FOUND: ${Python_FOUND}")
+message("Python_VERSION: ${Python_VERSION}")
+message("Python_Development_FOUND: ${Python_Development_FOUND}")
+message("Python_LIBRARIES: ${Python_LIBRARIES}")
+message("Python_INCLUDE_DIRS: ${Python_INCLUDE_DIRS}")
+
+# Ensure that numpy is installed and read its include dir
+exec_program(${Python_EXECUTABLE}
+        ARGS "-c \"import numpy; print(numpy.get_include())\""
+        OUTPUT_VARIABLE NUMPY_INCLUDE_DIR
+        RETURN_VALUE NUMPY_NOT_FOUND
+        )
+if (NUMPY_NOT_FOUND)
+    message(FATAL_ERROR "NumPy headers not found")
+endif ()
+
+# Ensure that pybind11 is installed and read its include dir
+exec_program(${Python_EXECUTABLE}
+        ARGS "-c \"import pybind11; print(pybind11.get_include())\""
+        OUTPUT_VARIABLE PYBIND11_INCLUDE_DIR
+        RETURN_VALUE PYBIND11_NOT_FOUND
+        )
+if (PYBIND11_NOT_FOUND)
+    message(FATAL_ERROR "PyBind11 headers not found")
+endif ()
+
+# Determine correct extension suffix
+exec_program(${Python_EXECUTABLE}
+        ARGS "-c \"import distutils.sysconfig; print(distutils.sysconfig.get_config_var('EXT_SUFFIX'))\""
+        OUTPUT_VARIABLE EXT_SUFFIX
+        RETURN_VALUE SUFFIX_NOT_FOUND
+        )
+if (SUFFIX_NOT_FOUND)
+    message(FATAL_ERROR "Extension suffix not found")
+endif ()
+
+#include_directories(src)
+add_subdirectory(src)
+
+#
+# Copy all python scripts to the build directory.
+#
+set(Python_SCRIPTS scattnlay/__init__.py scattnlay/main.py)
+
+foreach (_script ${Python_SCRIPTS})
+    configure_file(
+            ${PROJECT_SOURCE_DIR}/${_script}
+            ${PROJECT_BINARY_DIR}/${_script}
+            COPYONLY
+    )
+endforeach ()

Makefile

@@ -37,13 +37,13 @@ deb: source
 	# build the package
 	dpkg-buildpackage -i -I -rfakeroot
 
-standalone: scattnlay fieldnlay scattnlay-mp fieldnlay-mp
+standalone: scattnlay-dp fieldnlay-dp scattnlay-mp fieldnlay-mp
 
 # standalone programs with DP
-scattnlay: $(SRCDIR)/farfield.cc $(SRCDIR)/nmie.cc $(CXX_NMIE_HEADERS)
+scattnlay-dp: $(SRCDIR)/farfield.cc $(SRCDIR)/nmie.cc $(CXX_NMIE_HEADERS)
 	$(CXX) -DNDEBUG -O2 -Wall -std=c++11 $(SRCDIR)/farfield.cc $(SRCDIR)/nmie.cc  -lm -o scattnlay-dp $(CXXFLAGS) $(LDFLAGS)
 
-fieldnlay: $(SRCDIR)/nearfield.cc $(SRCDIR)/nmie.cc $(CXX_NMIE_HEADERS)
+fieldnlay-dp: $(SRCDIR)/nearfield.cc $(SRCDIR)/nmie.cc $(CXX_NMIE_HEADERS)
 	$(CXX) -DNDEBUG -O2 -Wall -std=c++11 $(SRCDIR)/nearfield.cc $(SRCDIR)/nmie.cc  -lm -o fieldnlay-dp $(CXXFLAGS) $(LDFLAGS)
 
 # standalone programs with MP

README.md

@@ -76,11 +76,11 @@ Binary files for Ubuntu and derivative distributions can be found at
 To install it you must configure the repository:
 ``` bash
 sudo add-apt-repository ppa:ovidio/scattering
-sudo apt-get update
+sudo apt update
 ```
 and then you simply install the package:
 ``` bash
-sudo apt-get install python-scattnlay
+sudo apt install python-scattnlay
 ```
 You can also install it from PyPi via
 ```bash

debian/control

@@ -2,16 +2,17 @@ Source: python-scattnlay
 Section: science
 Priority: optional
 Maintainer: Ovidio Peña Rodríguez <ovidio@bytesfall.com>
-Build-Depends: debhelper (>=7.0.0), dh-python, cdbs (>= 0.4.49), libboost-all-dev (>= 1.58.0), python-numpy (>= 1.0), python-all-dev, python-numpy-dev
-X-Python-Version: >=2.7
+Build-Depends: debhelper (>= 9), dh-python, cdbs (>= 0.4.49), libboost-all-dev (>= 1.65.1), python3-numpy (>= 1.0), python3-all-dev, python3-numpy-dev
+X-Python3-Version: >= 3.2
 Standards-Version: 3.8.3
 
-Package: python-scattnlay
+Package: python3-scattnlay
 Architecture: any
 Homepage: http://scattering.sourceforge.net/
-XB-Python-Version: ${python:Versions}
-Depends: ${shlibs:Depends}, ${misc:Depends}, ${python:Depends}, libboost-math (>= 1.58.0), python-numpy (>= 1.0)
-Description: Calculation of the scattering of EM radiation by a multilayered sphere
+XB-Python3-Version: ${python3:Versions}
+Depends: ${shlibs:Depends}, ${misc:Depends}, ${python3:Depends}, libboost-math1.65.1 (>= 1.65.1), python3-numpy (>= 1.0)
+Description: Python scattnlay (Python 3)
+ Calculation of the scattering of EM radiation by a multilayered sphere
  The Python version of scattnlay, a computer implementation of the algorithm
  for the calculation of electromagnetic radiation scattering by a multilayered
  sphere developed by Yang. It has been shown that the program is effective,

debian/rules

@@ -1,7 +1,7 @@
 #!/usr/bin/make -f
 # -*- makefile -*-
 
-export CFLAGS='-std=c++11' '-Wc++0x-compat'
+export CFLAGS = '-std=c++11' '-Wc++0x-compat'
 
 include /usr/share/cdbs/1/rules/debhelper.mk
 include /usr/share/cdbs/1/class/python-distutils.mk

examples/coating-flow.py

@@ -56,7 +56,7 @@ if __name__ == '__main__':
     args = parser.parse_args()
 
     for dirname in args.dirnames:
-        print "Calculating spectra for data file(s) in dir '%s'..." % (dirname)
+        print("Calculating spectra for data file(s) in dir '%s'..." % (dirname))
 
         wl = args.wl # cm
         if (args.rad is None):
@@ -78,11 +78,11 @@ if __name__ == '__main__':
 
         Rt = Rs + tc # cm
 
-        print "Wl = %.2f, Rs = %.2f, tc = %.2f, Rt = %.2f" % (wl, Rs, tc, Rt)
+        print("Wl = %.2f, Rs = %.2f, tc = %.2f, Rt = %.2f" % (wl, Rs, tc, Rt))
 
         ms = 1.0 + 40.0j
         for i, fname in enumerate(files):
-            print "Calculating spectra for file '%s'..." % (fname)
+            print("Calculating spectra for file '%s'..." % (fname))
 
             basename = os.path.splitext(fname)[0]
 
@@ -99,12 +99,12 @@ if __name__ == '__main__':
 
             x = 2.0*np.pi*np.array(r, dtype = np.float64)/wl
             m = np.array([ms] + nvalues[:, 1].tolist(), dtype = np.complex128)
-            print(x,m)
+            #print(x,m)
 
             factor = 2.91*x[0]/x[-1]
-            print factor
+            #print(factor)
             comment='PEC-'+basename
-            WL_units=''
+            WL_units='cm'
             #flow_total = 39
             # flow_total = 23 #SV False
             flow_total = 24
@@ -115,28 +115,23 @@ if __name__ == '__main__':
             #crossplane='XY'
 
             # Options to plot: Eabs, Habs, Pabs, angleEx, angleHy
-            #field_to_plot='Pabs'
+            field_to_plot='Pabs'
             #field_to_plot='Eabs'
             
-            field_to_plot='angleEx'
+            #field_to_plot='angleEx'
             #field_to_plot='angleHy'
-            print "x =", x
-            print "m =", m
+            #print("x =", x)
+            #print("m =", m)
 
             import matplotlib.pyplot as plt
             plt.rcParams.update({'font.size': 16})
             fig, axs = plt.subplots(1,1)#, sharey=True, sharex=True)
             fig.tight_layout()
             fieldplot(fig, axs, x,m, wl, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total,
-                      subplot_label=' ',is_flow_extend=False
-                      , outline_width=1.5
-                      , pl=0 #PEC layer starts the design
+                      subplot_label=' ', outline_width=1.5,
+                      pl=0 #PEC layer starts the design
                       )
-            # fieldplot(fig, axs, x[0],m[0], wl, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total,
-            #           subplot_label=' ' ,is_flow_extend=False
-            #           , outline_width=1.5
-            #           , pl=0 #PEC layer starts the design
-            #           )
+
             fig.subplots_adjust(hspace=0.3, wspace=-0.1)
             plt.savefig(comment+"-R"+str(int(round(x[-1]*wl/2.0/np.pi)))+"-"+crossplane+"-"
                         +field_to_plot+".pdf",pad_inches=0.02, bbox_inches='tight')
@@ -145,5 +140,5 @@ if __name__ == '__main__':
             plt.close()
 
 
-        print "Done!!"
+        print("Done!!")
 

examples/coating-flow.sh

@@ -5,4 +5,4 @@
 #python coating-flow.py -w 3.75 -t 0.62 -f index-sv.dat -n 1501
 # python coating-flow.py -w 3.75 -t 2.40 -f index-ch.dat -n 501
 
-python coating-flow.py -w 0.532 -t 0.128 -f index-in-glass.dat -n 151
+python3 coating-flow.py -w 0.532 -t 0.128 -f index-in-glass.dat -n 151

examples/field-Ag-flow.py

@@ -76,7 +76,7 @@ comment = 'bulk-WL' + str(WL) + 'nm_r' + str(core_r) + 'nm_epsilon' + str(epsilo
 WL_units = 'nm'
 npts = 251
 factor = 2.1
-flow_total = 9
+flow_total = 41
 # flow_total = 21
 # flow_total = 0
 crossplane = 'XZ'
@@ -92,8 +92,9 @@ import matplotlib.pyplot as plt
 
 fig, axs = plt.subplots(1, 1)  # , sharey=True, sharex=True)
 fig.tight_layout()
-fieldplot(fig, axs, x, m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total,
-          subplot_label=' ', is_flow_extend=False)
+fieldplot(fig, axs, x, m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor,
+          flow_total, density=50.0, maxlength=40.0, arrowstyle='-',
+          subplot_label=' ', draw_shell=True)
 
 # fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total, is_flow_extend=False)
 

examples/field-SiAgSi-flow.py

@@ -91,9 +91,10 @@ factor=2.2
 flow_total = 21
 
 
-crossplane='XZ'
+#crossplane='XZ'
 #crossplane='YZ'
 #crossplane='XY'
+crossplane='XYZ'
 
 # Options to plot: Eabs, Habs, Pabs, angleEx, angleHy
 field_to_plot='Eabs'
@@ -105,9 +106,9 @@ import matplotlib.pyplot as plt
 fig, axs = plt.subplots(1,1)#, sharey=True, sharex=True)
 fig.tight_layout()
 fieldplot(fig, axs, x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total,
-          subplot_label=' ',is_flow_extend=False, outline_width=1.5)
+          subplot_label=' ', outline_width=1.5, draw_shell=True)
 
-#fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total, is_flow_extend=False)
+#fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total)
 
 # for ax in axs:
 #     ax.locator_params(axis='x',nbins=5)

examples/fieldplot.py

@@ -28,153 +28,41 @@
 
 # Several functions to plot field and streamlines (power flow lines).
 
-import scattnlay
-from scattnlay import fieldnlay
-from scattnlay import scattnlay
+from scattnlay import fieldnlay, scattnlay
 import numpy as np
-import cmath
-
-
-def unit_vector(vector):
-    """ Returns the unit vector of the vector.  """
-    return vector / np.linalg.norm(vector)
-
-
-def angle_between(v1, v2):
-    """ Returns the angle in radians between vectors 'v1' and 'v2'::
-
-            >>> angle_between((1, 0, 0), (0, 1, 0))
-            1.5707963267948966
-            >>> angle_between((1, 0, 0), (1, 0, 0))
-            0.0
-            >>> angle_between((1, 0, 0), (-1, 0, 0))
-            3.141592653589793
-    """
-    v1_u = unit_vector(v1)
-    v2_u = unit_vector(v2)
-    angle = np.arccos(np.dot(v1_u, v2_u))
-    if np.isnan(angle):
-        if (v1_u == v2_u).all():
-            return 0.0
-        else:
-            return np.pi
-    return angle
-###############################################################################
-
-
-def GetFlow3D(x0, y0, z0, max_length, max_angle, x, m, pl):
-    # Initial position
-    flow_x = [x0]
-    flow_y = [y0]
-    flow_z = [z0]
-    max_step = x[-1] / 3
-    min_step = x[0] / 2000
-#    max_step = min_step
-    step = min_step * 2.0
-    if max_step < min_step:
-        max_step = min_step
-    terms, E, H = fieldnlay(np.array([x]), np.array([m]), np.array([flow_x[-1]]), np.array([flow_y[-1]]), np.array([flow_z[-1]]), pl=pl)
-    Ec, Hc = E[0, 0, :], H[0, 0, :]
-    S = np.cross(Ec, Hc.conjugate()).real
-    Snorm_prev = S / np.linalg.norm(S)
-    Sprev = S
-    length = 0
-    dpos = step
-    count = 0
-    while length < max_length:
-        count = count + 1
-        if (count > 4000):  # Limit length of the absorbed power streamlines
-            break
-        if step < max_step:
-            step = step * 2.0
-        r = np.sqrt(flow_x[-1]**2 + flow_y[-1]**2 + flow_z[-1]**2)
-        while step > min_step:
-            # Evaluate displacement from previous poynting vector
-            dpos = step
-            dx = dpos * Snorm_prev[0]
-            dy = dpos * Snorm_prev[1]
-            dz = dpos * Snorm_prev[2]
-            # Test the next position not to turn\chang size for more than
-            # max_angle
-            coord = np.vstack((np.array([flow_x[-1] + dx]), np.array([flow_y[-1] + dy]),
-                               np.array([flow_z[-1] + dz]))).transpose()
-            terms, E, H = fieldnlay(np.array([x]), np.array([m]), np.array([flow_x[-1] + dx]),
-                                    np.array([flow_y[-1] + dy]), np.array([flow_z[-1] + dz]), pl=pl)
-            Ec, Hc = E[0, 0, :], H[0, 0, :]
-            Eth = max(np.absolute(Ec)) / 1e10
-            Hth = max(np.absolute(Hc)) / 1e10
-            for i in range(0, len(Ec)):
-                if abs(Ec[i]) < Eth:
-                    Ec[i] = 0 + 0j
-                if abs(Hc[i]) < Hth:
-                    Hc[i] = 0 + 0j
-            S = np.cross(Ec, Hc.conjugate()).real
-            if not np.isfinite(S).all():
-                break
-            Snorm = S / np.linalg.norm(S)
-            diff = (S - Sprev) / max(np.linalg.norm(S), np.linalg.norm(Sprev))
-            if np.linalg.norm(diff) < max_angle:
-                # angle = angle_between(Snorm, Snorm_prev)
-                # if abs(angle) < max_angle:
-                break
-            step = step / 2.0
-        # 3. Save result
-        Sprev = S
-        Snorm_prev = Snorm
-        dx = dpos * Snorm_prev[0]
-        dy = dpos * Snorm_prev[1]
-        dz = dpos * Snorm_prev[2]
-        length = length + step
-        flow_x.append(flow_x[-1] + dx)
-        flow_y.append(flow_y[-1] + dy)
-        flow_z.append(flow_z[-1] + dz)
-    return np.array(flow_x), np.array(flow_y), np.array(flow_z)
 
 
 ###############################################################################
-def GetField(crossplane, npts, factor, x, m, pl):
+def GetCoords(crossplane, npts, factor, x):
     """
     crossplane: XZ, YZ, XY, or XYZ (half is XZ, half is YZ)
     npts: number of point in each direction
     factor: ratio of plotting size to outer size of the particle
     x: size parameters for particle layers
-    m: relative index values for particle layers
     """
     scan = np.linspace(-factor*x[-1], factor*x[-1], npts)
     zero = np.zeros(npts*npts, dtype = np.float64)
 
     if crossplane=='XZ':
         coordX, coordZ = np.meshgrid(scan, scan)
-        coordX.resize(npts * npts)
-        coordZ.resize(npts * npts)
+        coordX.resize(npts*npts)
+        coordZ.resize(npts*npts)
         coordY = zero
-        coordPlot1 = coordX
-        coordPlot2 = coordZ
     elif crossplane == 'YZ':
         coordY, coordZ = np.meshgrid(scan, scan)
-        coordY.resize(npts * npts)
-        coordZ.resize(npts * npts)
+        coordY.resize(npts*npts)
+        coordZ.resize(npts*npts)
         coordX = zero
-        coordPlot1 = coordY
-        coordPlot2 = coordZ
     elif crossplane == 'XY':
         coordX, coordY = np.meshgrid(scan, scan)
-        coordX.resize(npts * npts)
-        coordY.resize(npts * npts)
+        coordX.resize(npts*npts)
+        coordY.resize(npts*npts)
         coordZ = zero
-        coordPlot1 = coordY
-        coordPlot2 = coordX
-    elif crossplane=='XYZ':
+    elif crossplane=='XYZ': # Upper half: XZ; Lower half: YZ
         coordX, coordZ = np.meshgrid(scan, scan)
         coordY, coordZ = np.meshgrid(scan, scan)
-        coordPlot1, coordPlot2 = np.meshgrid(scan, scan)
-        coordPlot1.resize(npts * npts)
-        coordPlot2.resize(npts * npts)
-        half=npts//2
-        # coordX = np.copy(coordX)
-        # coordY = np.copy(coordY)
-        coordX[:,:half]=0
-        coordY[:,half:]=0
+        coordX[:, scan<0] = 0
+        coordY[:, scan>=0] = 0
         coordX.resize(npts*npts)
         coordY.resize(npts*npts)
         coordZ.resize(npts*npts)
@@ -183,88 +71,102 @@ def GetField(crossplane, npts, factor, x, m, pl):
         import sys
         sys.exit()
 
-    terms, E, H = fieldnlay(np.array([x]), np.array([m]), coordX, coordY, coordZ, pl=pl)
-    Ec = E[0, :, :]
-    Hc = H[0, :, :]
-    P = np.array(list(map(lambda n: np.linalg.norm(np.cross(Ec[n],
-                                                            np.conjugate(Hc[n])
-                                                            # Hc[n]
-                                                            )).real,
-                     range(0, len(E[0])))))
-    print(P)
-    # for n in range(0, len(E[0])):
-    #     P.append(np.linalg.norm( np.cross(Ec[n], np.conjugate(Hc[n]) ).real/2 ))
-    return Ec, Hc, P, coordPlot1, coordPlot2
+    return coordX, coordY, coordZ, scan
+
+
+###############################################################################
+def GetField(crossplane, npts, factor, x, m, pl):
+    """
+    crossplane: XZ, YZ, XY, or XYZ (half is XZ, half is YZ)
+    npts: number of point in each direction
+    factor: ratio of plotting size to outer size of the particle
+    x: size parameters for particle layers
+    m: relative index values for particle layers
+    """
+    coordX, coordY, coordZ, scan = GetCoords(crossplane, npts, factor, x)
+
+    terms, E, H = fieldnlay(x, m, coordX, coordY, coordZ, pl=pl)
+    if len(E.shape) > 2:
+        E = E[0, :, :]
+        H = H[0, :, :]
+
+    S = np.cross(E, np.conjugate(H)).real
+    print(S)
+
+    if crossplane=='XZ':
+        Sx = np.resize(S[:, 2], (npts, npts)).T
+        Sy = np.resize(S[:, 0], (npts, npts)).T
+    elif crossplane == 'YZ':
+        Sx = np.resize(S[:, 2], (npts, npts)).T
+        Sy = np.resize(S[:, 1], (npts, npts)).T
+    elif crossplane == 'XY':
+        Sx = np.resize(S[:, 1], (npts, npts)).T
+        Sy = np.resize(S[:, 0], (npts, npts)).T
+    elif crossplane=='XYZ': # Upper half: XZ; Lower half: YZ
+        Sx = np.resize(S[:, 2], (npts, npts)).T
+        Sy = np.resize(S[:, 0], (npts, npts)).T
+        Sy[scan<0] = np.resize(S[:, 1], (npts, npts)).T[scan<0]
+    else:
+        print("Unknown crossplane")
+        import sys
+        sys.exit()
+
+    return E, H, S, scan, Sx, Sy
 ###############################################################################
 
 
 def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
-              field_to_plot='Pabs', npts=101, factor=2.1, flow_total=11,
-              is_flow_extend=True, pl=-1, outline_width=1, subplot_label=' '):
-    # print(fig, ax, x, m, WL, comment, WL_units, crossplane,
-    #       field_to_plot, npts, factor, flow_total,
-    #       is_flow_extend, pl, outline_width, subplot_label)
-    Ec, Hc, P, coordX, coordZ = GetField(crossplane, npts, factor, x, m, pl)
-    Er = np.absolute(Ec)
-    Hr = np.absolute(Hc)
+              field_to_plot='Pabs', npts=101, factor=2.1,
+              flow_total=11, density=20, minlength=0.1, maxlength=4.0,
+              arrowstyle='-|>', arrowsize=1.0,
+              pl=-1, draw_shell=False, outline_width=1, subplot_label=' '):
+
+    E, H, S, scan, Sx, Sy = GetField(crossplane, npts, factor, x, m, pl)
+    Er = np.absolute(E)
+    Hr = np.absolute(H)
     try:
         from matplotlib import cm
         from matplotlib.colors import LogNorm
 
         if field_to_plot == 'Pabs':
-            Eabs_data = np.resize(P, (npts, npts)).T
             label = r'$\operatorname{Re}(E \times H^*)$'
+            data = np.resize(np.linalg.norm(np.cross(E, np.conjugate(H)), axis=1).real, (npts, npts)).T
         elif field_to_plot == 'Eabs':
-            Eabs = np.sqrt(Er[:, 0]**2 + Er[:, 1]**2 + Er[:, 2]**2)
             label = r'$|E|$'
-            # Eabs = np.real(Hc[:, 0])
-            # label = r'$Re(H_x)$'
-            # Eabs = np.real(Hc[:, 1])
-            # label = r'$Re(H_y)$'
-            # Eabs = np.real(Ec[:, 1])
-            # label = r'$Re(E_y)$'
-            # Eabs = np.real(Ec[:, 0])
-            # label = r'$Re(E_x)$'
-            Eabs_data = np.resize(Eabs, (npts, npts)).T
+            Eabs = np.sqrt(Er[:, 0]**2 + Er[:, 1]**2 + Er[:, 2]**2)
+            data = np.resize(Eabs, (npts, npts)).T
         elif field_to_plot == 'Habs':
-            Habs = np.sqrt(Hr[:, 0]**2 + Hr[:, 1]**2 + Hr[:, 2]**2)
-            Habs = 376.730313667 * Habs # scale to free space impedance
-            Eabs_data = np.resize(Habs, (npts, npts)).T
             label = r'$|H|$'
+            Habs = np.sqrt(Hr[:, 0]**2 + Hr[:, 1]**2 + Hr[:, 2]**2)
+            Habs = 376.730313667*Habs # scale to free space impedance
+            data = np.resize(Habs, (npts, npts)).T
         elif field_to_plot == 'angleEx':
-            Eangle = np.angle(Ec[:, 0]) / np.pi * 180
-            Eabs_data = np.resize(Eangle, (npts, npts)).T
             label = r'$arg(E_x)$'
+            Eangle = np.angle(E[:, 0])/np.pi*180
+            data = np.resize(Eangle, (npts, npts)).T
         elif field_to_plot == 'angleHy':
-            Hangle = np.angle(Hc[:, 1]) / np.pi * 180
-            Eabs_data = np.resize(Hangle, (npts, npts)).T
             label = r'$arg(H_y)$'
+            Hangle = np.angle(H[:, 1])/np.pi*180
+            data = np.resize(Hangle, (npts, npts)).T
 
         # Rescale to better show the axes
-        scale_x = np.linspace(
-            min(coordX) * WL / 2.0 / np.pi, max(coordX) * WL / 2.0 / np.pi, npts)
-        scale_z = np.linspace(
-            min(coordZ) * WL / 2.0 / np.pi, max(coordZ) * WL / 2.0 / np.pi, npts)
+        scale = scan*WL/2.0/np.pi
 
         # Define scale ticks
-        min_tick = np.amin(Eabs_data[~np.isnan(Eabs_data)])
-        #min_tick = 0.1
-        max_tick = np.amax(Eabs_data[~np.isnan(Eabs_data)])
-        #max_tick = 60
+        min_tick = np.amin(data[~np.isnan(data)])
+        max_tick = np.amax(data[~np.isnan(data)])
+
         scale_ticks = np.linspace(min_tick, max_tick, 5)
-        #scale_ticks = np.power(10.0, np.linspace(np.log10(min_tick), np.log10(max_tick), 6))
-        #scale_ticks = [0.1,0.3,1,3,10, max_tick]
-        # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
+
         ax.set_title(label)
         # build a rectangle in axes coords
         ax.annotate(subplot_label, xy=(0.0, 1.1), xycoords='axes fraction',  # fontsize=10,
                     horizontalalignment='left', verticalalignment='top')
-        # ax.text(right, top, subplot_label,
-        #         horizontalalignment='right',
-        #         verticalalignment='bottom',
-        #         transform=ax.transAxes)
-        cax = ax.imshow(Eabs_data, interpolation='nearest', cmap=cm.jet,
-                        origin='lower', vmin=min_tick, vmax=max_tick, extent=(min(scale_x), max(scale_x), min(scale_z), max(scale_z))
+
+        # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
+        cax = ax.imshow(data, interpolation='nearest', cmap=cm.jet,
+                        origin='lower', vmin=min_tick, vmax=max_tick,
+                        extent=(min(scale), max(scale), min(scale), max(scale))
                         # ,norm = LogNorm()
                         )
         ax.axis("image")
@@ -276,108 +178,47 @@ def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
             cbar.ax.set_yticklabels(['%3.0f' % (a) for a in scale_ticks])
         else:
             cbar.ax.set_yticklabels(['%g' % (a) for a in scale_ticks])
-        # pos = list(cbar.ax.get_position().bounds)
-        #fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
-        lp2 = -10.0
-        lp1 = -1.0
+
         if crossplane == 'XZ':
-            ax.set_xlabel('Z, ' + WL_units, labelpad=lp1)
-            ax.set_ylabel('X, ' + WL_units, labelpad=lp2)
+            ax.set_xlabel('Z (%s)' % (WL_units))
+            ax.set_ylabel('X (%s)' % (WL_units))
         elif crossplane == 'YZ':
-            ax.set_xlabel('Z, ' + WL_units, labelpad=lp1)
-            ax.set_ylabel('Y, ' + WL_units, labelpad=lp2)
+            ax.set_xlabel('Z (%s)' % (WL_units))
+            ax.set_ylabel('Y (%s)' % (WL_units))
         elif crossplane=='XYZ':
-            ax.set_xlabel(r'$Z,\lambda$'+WL_units)
-            ax.set_ylabel(r'$Y:X,\lambda$'+WL_units)
+            ax.set_xlabel('Z (%s)' % (WL_units))
+            ax.set_ylabel('Y(<0):X(>0) (%s)' % (WL_units))
+
+            # draw a line to separate both planes
+            ax.axhline(linewidth=outline_width, color='black')
         elif crossplane == 'XY':
-            ax.set_xlabel('X, ' + WL_units, labelpad=lp1)
-            ax.set_ylabel('Y, ' + WL_units, labelpad=lp2)
-        # # This part draws the nanoshell
-        from matplotlib import patches
-        from matplotlib.path import Path
-        for xx in x:
-            r = xx * WL / 2.0 / np.pi
-            s1 = patches.Arc((0, 0), 2.0 * r, 2.0 * r,  angle=0.0, zorder=1.8,
-                             theta1=0.0, theta2=360.0, linewidth=outline_width, color='black')
-            ax.add_patch(s1)
-        #
-        # for flow in range(0,flow_total):
-        #     flow_x, flow_z = GetFlow(scale_x, scale_z, Ec, Hc,
-        #                              min(scale_x)+flow*(scale_x[-1]-scale_x[0])/(flow_total-1),
-        #                              min(scale_z),
-        #                              #0.0,
-        #                              npts*16)
-        #     verts = np.vstack((flow_z, flow_x)).transpose().tolist()
-        #     #codes = [Path.CURVE4]*len(verts)
-        #     codes = [Path.LINETO]*len(verts)
-        #     codes[0] = Path.MOVETO
-        #     path = Path(verts, codes)
-        #     patch = patches.PathPatch(path, facecolor='none', lw=1, edgecolor='yellow')
-        #     ax.add_patch(patch)
-        if (not crossplane == 'XY') and flow_total > 0:
+            ax.set_xlabel('X (%s)' % (WL_units))
+            ax.set_ylabel('Y (%s)' % (WL_units))
 
+        if draw_shell:
+            # Draw the nanoshell
+            from matplotlib import patches
             from matplotlib.path import Path
-            scanSP = np.linspace(-factor * x[-1], factor * x[-1], npts)
-            min_SP = -factor * x[-1]
-            step_SP = 2.0 * factor * x[-1] / (flow_total - 1)
-            x0, y0, z0 = 0, 0, 0
-            max_length = factor * x[-1] * 10
-            # max_length=factor*x[-1]*5
-            max_angle = np.pi / 160
-            if is_flow_extend:
-                rg = range(0, flow_total * 5 + 1)
-            else:
-                rg = range(0, flow_total)
-            for flow in rg:
-                if is_flow_extend:
-                    f = min_SP*2 + flow*step_SP
-                else:
-                    f = min_SP + flow*step_SP
-                if crossplane=='XZ':
-                    x0 = f
-                elif crossplane=='YZ':
-                    y0 = f
-                elif crossplane=='XYZ':
-                    x0 = 0
-                    y0 = 0
-                    if f > 0:
-                        x0 = f
-                    else:
-                        y0 = f
-                z0 = min_SP
-                    # x0 = x[-1]/20
-                flow_xSP, flow_ySP, flow_zSP = GetFlow3D(
-                    x0, y0, z0, max_length, max_angle, x, m, pl)
-                if crossplane == 'XZ':
-                    flow_z_plot = flow_zSP * WL / 2.0 / np.pi
-                    flow_f_plot = flow_xSP * WL / 2.0 / np.pi
-                elif crossplane == 'YZ':
-                    flow_z_plot = flow_zSP * WL / 2.0 / np.pi
-                    flow_f_plot = flow_ySP * WL / 2.0 / np.pi
-                elif crossplane=='XYZ':
-                    if f > 0:
-                        flow_z_plot = flow_zSP*WL/2.0/np.pi
-                        flow_f_plot = flow_xSP*WL/2.0/np.pi
-                    else:
-                        flow_z_plot = flow_zSP*WL/2.0/np.pi
-                        flow_f_plot = flow_ySP*WL/2.0/np.pi
-
-                verts = np.vstack(
-                    (flow_z_plot, flow_f_plot)).transpose().tolist()
-                codes = [Path.LINETO] * len(verts)
-                codes[0] = Path.MOVETO
-                path = Path(verts, codes)
-                #patch = patches.PathPatch(path, facecolor='none', lw=0.2, edgecolor='white',zorder = 2.7)
-                patch = patches.PathPatch(
-                    path, facecolor='none', lw=outline_width, edgecolor='white', zorder=1.9, alpha=0.7)
-                # patch = patches.PathPatch(
-                #     path, facecolor='none', lw=0.7, edgecolor='white', zorder=1.9, alpha=0.7)
-                ax.add_patch(patch)
-#                ax.plot(flow_z_plot, flow_f_plot, 'x', ms=2, mew=0.1,
-#                        linewidth=0.5, color='k', fillstyle='none')
+            for xx in x:
+                r = xx*WL/2.0/np.pi
+                s1 = patches.Arc((0, 0), 2.0*r, 2.0*r,  angle=0.0, zorder=1.8,
+                                 theta1=0.0, theta2=360.0, linewidth=outline_width, color='black')
+                ax.add_patch(s1)
 
+        # Draw flow lines
+        if (not crossplane == 'XY') and flow_total > 0:
+            margin = 0.98
+            points = np.vstack((margin*scale.min()*np.ones(flow_total),
+                                np.linspace(margin*scale.min(),
+                                            margin*scale.max(), flow_total))).transpose()
+
+            # Plot the streamlines with an appropriate colormap and arrow style
+            ax.streamplot(scale, scale, Sx, Sy,
+                          start_points=points, integration_direction='both',
+                          density=density, minlength=minlength, maxlength=maxlength,
+                          linewidth=outline_width, color='white',
+                          arrowstyle=arrowstyle, arrowsize=arrowsize)
     finally:
-        terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(
-            np.array([x]), np.array([m]))
+        terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
         print("Qsca = " + str(Qsca))
     #

requirements.txt

@@ -0,0 +1,6 @@
+scipy~=1.3.3
+numpy~=1.17.4
+matplotlib~=3.1.2
+PyYAML~=5.3.1
+pybind11~=2.4.3
+setuptools~=45.2.0

scattnlay/__init__.py

@@ -30,4 +30,4 @@
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-from scattnlay.main import scattcoeffs, scattnlay, fieldnlay
+from scattnlay.main import scattcoeffs, expancoeffs, scattnlay, fieldnlay

scattnlay/main.py

@@ -68,6 +68,10 @@ def scattcoeffs(x, m, nmax=-1, pl=-1, mp=False):
     elif len(x.shape) != 2:
         raise ValueError('The size parameter (x) should be a 1-D or 2-D NumPy array.')
 
+    # Repeat the same m for all wavelengths
+    if len(m.shape) == 1:
+        m = np.repeat(m[np.newaxis, :], x.shape[0], axis=0)
+
     if nmax == -1:
         nstore = 0
     else:
@@ -78,12 +82,7 @@ def scattcoeffs(x, m, nmax=-1, pl=-1, mp=False):
     bn = np.zeros((0, nstore), dtype=complex)
 
     for i, xi in enumerate(x):
-        if len(m.shape) == 1:
-            mi = m
-        else:
-            mi = m[i]
-
-        terms[i], a, b = scattcoeffs_(xi, mi, nmax=nmax, pl=pl)
+        terms[i], a, b = scattcoeffs_(xi, m[i], nmax=nmax, pl=pl)
 
         if terms[i] > nstore:
             nstore = terms[i]
@@ -96,6 +95,75 @@ def scattcoeffs(x, m, nmax=-1, pl=-1, mp=False):
     return terms, an, bn
 #scattcoeffs()
 
+def expancoeffs(x, m, nmax=-1, pl=-1, mp=False):
+    """
+    expancoeffs(x, m[, nmax, pl, mp])
+
+    Calculate the scattering coefficients required to calculate both the
+    near- and far-field parameters.
+
+        x: Size parameters (1D or 2D ndarray)
+        m: Relative refractive indices (1D or 2D ndarray)
+        nmax: Maximum number of multipolar expansion terms to be used for the
+              calculations. Only use it if you know what you are doing, otherwise
+              set this parameter to -1 and the function will calculate it.
+        pl: Index of PEC layer. If there is none just send -1.
+        mp: Use multiple (True) or double (False) precision.
+
+    Returns: (terms, an, bn, cn, dn)
+    with
+        terms: Number of multipolar expansion terms used for the calculations
+        an, bn, cn, dn: Complex expansion coefficients of each layer
+    """
+
+    if mp:
+        from scattnlay_mp import expancoeffs as expancoeffs_
+    else:
+        from scattnlay_dp import expancoeffs as expancoeffs_
+
+    if len(m.shape) != 1 and len(m.shape) != 2:
+        raise ValueError('The relative refractive index (m) should be a 1-D or 2-D NumPy array.')
+    if len(x.shape) == 1:
+        if len(m.shape) == 1:
+            return expancoeffs_(x, m, nmax=nmax, pl=pl)
+        else:
+            raise ValueError('The number of of dimensions for the relative refractive index (m) and for the size parameter (x) must be equal.')
+    elif len(x.shape) != 2:
+        raise ValueError('The size parameter (x) should be a 1-D or 2-D NumPy array.')
+
+    # Repeat the same m for all wavelengths
+    if len(m.shape) == 1:
+        m = np.repeat(m[np.newaxis, :], x.shape[0], axis=0)
+
+    if nmax == -1:
+        nstore = 0
+    else:
+        nstore = nmax
+
+    terms = np.zeros((x.shape[0]), dtype=int)
+    an = np.zeros((0, x.shape[1], nstore), dtype=complex)
+    bn = np.zeros((0, x.shape[1], nstore), dtype=complex)
+    cn = np.zeros((0, x.shape[1], nstore), dtype=complex)
+    dn = np.zeros((0, x.shape[1], nstore), dtype=complex)
+
+    for i, xi in enumerate(x):
+        terms[i], a, b, c, d = expancoeffs_(xi, m[i], nmax=nmax, pl=pl)
+
+        if terms[i] > nstore:
+            nstore = terms[i]
+            an.resize((an.shape[0], an.shape[1], nstore))
+            bn.resize((bn.shape[0], bn.shape[1], nstore))
+            cn.resize((cn.shape[0], cn.shape[1], nstore))
+            dn.resize((dn.shape[0], dn.shape[1], nstore))
+
+        an = np.vstack((an, a))
+        bn = np.vstack((bn, b))
+        cn = np.vstack((cn, c))
+        dn = np.vstack((dn, d))
+
+    return terms, an, bn, cn, dn
+#expancoeffs()
+
 
 def scattnlay(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
     """
@@ -142,6 +210,10 @@ def scattnlay(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
     if len(theta.shape) != 1:
         raise ValueError('The scattering angles (theta) should be a 1-D NumPy array.')
 
+    # Repeat the same m for all wavelengths
+    if len(m.shape) == 1:
+        m = np.repeat(m[np.newaxis, :], x.shape[0], axis=0)
+
     terms = np.zeros((x.shape[0]), dtype=int)
     Qext = np.zeros((x.shape[0]), dtype=float)
     Qsca = np.zeros((x.shape[0]), dtype=float)
@@ -154,12 +226,7 @@ def scattnlay(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
     S2 = np.zeros((x.shape[0], theta.shape[0]), dtype=complex)
 
     for i, xi in enumerate(x):
-        if len(m.shape) == 1:
-            mi = m
-        else:
-            mi = m[i]
-
-        terms[i], Qext[i], Qsca[i], Qabs[i], Qbk[i], Qpr[i], g[i], Albedo[i], S1[i], S2[i] = scattnlay_(xi, mi, theta, nmax=nmax, pl=pl)
+        terms[i], Qext[i], Qsca[i], Qabs[i], Qbk[i], Qpr[i], g[i], Albedo[i], S1[i], S2[i] = scattnlay_(xi, m[i], theta, nmax=nmax, pl=pl)
 
     return terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2
 #scattnlay()
@@ -174,11 +241,11 @@ def fieldnlay(x, m, xp, yp, zp, nmax=-1, pl=-1, mp=False):
         x: Size parameters (1D or 2D ndarray)
         m: Relative refractive indices (1D or 2D ndarray)
         xp: Array containing all X coordinates to calculate the complex
-            electric and magnetic fields (1D ndarray)
+            electric and magnetic fields (1D* ndarray)
         yp: Array containing all Y coordinates to calculate the complex
-            electric and magnetic fields (1D ndarray)
+            electric and magnetic fields (1D* ndarray)
         zp: Array containing all Z coordinates to calculate the complex
-            electric and magnetic fields (1D ndarray)
+            electric and magnetic fields (1D* ndarray)
         nmax: Maximum number of multipolar expansion terms to be used for the
               calculations. Only use it if you know what you are doing.
         pl: Index of PEC layer. If there is none just send -1.
@@ -188,6 +255,9 @@ def fieldnlay(x, m, xp, yp, zp, nmax=-1, pl=-1, mp=False):
     with
         terms: Number of multipolar expansion terms used for the calculations
         E, H: Complex electric and magnetic field at the provided coordinates
+
+    *Note: We assume that the coordinates are referred to the first wavelength
+           (or structure) and correct it for the following ones
     """
 
     if mp:
@@ -205,17 +275,18 @@ def fieldnlay(x, m, xp, yp, zp, nmax=-1, pl=-1, mp=False):
     elif len(x.shape) != 2:
         raise ValueError('The size parameter (x) should be a 1-D or 2-D NumPy array.')
 
+    # Repeat the same m for all wavelengths
+    if len(m.shape) == 1:
+        m = np.repeat(m[np.newaxis, :], x.shape[0], axis=0)
+
     terms = np.zeros((x.shape[0]), dtype=int)
     E = np.zeros((x.shape[0], xp.shape[0], 3), dtype=complex)
     H = np.zeros((x.shape[0], xp.shape[0], 3), dtype=complex)
 
     for i, xi in enumerate(x):
-        if len(m.shape) == 1:
-            mi = m
-        else:
-            mi = m[i]
-
-        terms[i], E[i], H[i] = fieldnlay_(xi, mi, xp, yp, zp, nmax=nmax, pl=pl)
+        # (2020/05/12) We assume that the coordinates are referred to the first wavelength
+        #              (or structure) and correct it for the following ones
+        terms[i], E[i], H[i] = fieldnlay_(xi, m[i], xp*xi[-1]/x[0, -1], yp*xi[-1]/x[0, -1], zp*xi[-1]/x[0, -1], nmax=nmax, pl=pl)
 
     return terms, E, H
 #fieldnlay()

setup.py

@@ -43,33 +43,32 @@ from setuptools.extension import Extension
 import numpy as np
 import pybind11 as pb
 
-setup(name = __mod__,
-      version = __version__,
-      description = __title__,
-      long_description="""The Python version of scattnlay, a computer implementation of the algorithm for the calculation of electromagnetic \
-radiation scattering by a multilayered sphere developed by Yang. It has been shown that the program is effective, \
-resulting in very accurate values of scattering efficiencies for a wide range of size parameters, which is a \
-considerable improvement over previous implementations of similar algorithms. For details see: \
-O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354.""",
-      author = __author__,
-      author_email = __email__,
-      maintainer = __author__,
-      maintainer_email = __email__,
-      keywords = ['Mie scattering', 'Multilayered sphere', 'Efficiency factors', 'Cross-sections'],
-      url = __url__,
-      download_url = __download_url__,
-      license = 'GPL',
-      platforms = 'any',
-      packages = ['scattnlay'],#, 'scattnlay_dp', 'scattnlay_mp'],
-      ext_modules = [Extension("scattnlay_dp",
-                               ["src/nmie.cc", "src/nmie-pybind11.cc", "src/pb11_wrapper.cc"],
-                               language = "c++",
-                               include_dirs = [np.get_include(), pb.get_include()],
-                               extra_compile_args=['-std=c++11']),
-                     Extension("scattnlay_mp",
-                               ["src/nmie.cc", "src/nmie-pybind11.cc", "src/pb11_wrapper.cc"],
-                               language = "c++",
-                               include_dirs = [np.get_include(), pb.get_include()],
-                               extra_compile_args=['-std=c++11', '-DMULTI_PRECISION=100'])]
-)
-
+setup(name=__mod__,
+      version=__version__,
+      description=__title__,
+      long_description="""The Python version of scattnlay, a computer implementation of the algorithm for the \
+calculation of electromagnetic radiation scattering by a multilayered sphere developed by Yang. It has been \
+shown that the program is effective, resulting in very accurate values of scattering efficiencies for a wide \
+range of size parameters, which is a considerable improvement over previous implementations of similar algorithms. \
+For details see: O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354.""",
+      author=__author__,
+      author_email=__email__,
+      maintainer=__author__,
+      maintainer_email=__email__,
+      keywords=['Mie scattering', 'Multilayered sphere', 'Efficiency factors', 'Cross-sections'],
+      url=__url__,
+      download_url=__download_url__,
+      license='GPL',
+      platforms='any',
+      packages=['scattnlay'],  # , 'scattnlay_dp', 'scattnlay_mp'],
+      ext_modules=[Extension("scattnlay_dp",
+                             ["src/nmie.cc", "src/nmie-pybind11.cc", "src/pb11_wrapper.cc"],
+                             language="c++",
+                             include_dirs=[np.get_include(), pb.get_include()],
+                             extra_compile_args=['-std=c++11']),
+                   Extension("scattnlay_mp",
+                             ["src/nmie.cc", "src/nmie-pybind11.cc", "src/pb11_wrapper.cc"],
+                             language="c++",
+                             include_dirs=[np.get_include(), pb.get_include()],
+                             extra_compile_args=['-std=c++11', '-DMULTI_PRECISION=100'])]
+      )

src/CMakeLists.txt

@@ -0,0 +1,57 @@
+# Sources for python extension
+set(_scattnlay_python_sources
+        ${CMAKE_CURRENT_LIST_DIR}/nmie.hpp
+        ${CMAKE_CURRENT_LIST_DIR}/nmie.cc
+        ${CMAKE_CURRENT_LIST_DIR}/nmie-pybind11.hpp
+        ${CMAKE_CURRENT_LIST_DIR}/nmie-pybind11.cc
+        ${CMAKE_CURRENT_LIST_DIR}/nmie-precision.hpp
+        ${CMAKE_CURRENT_LIST_DIR}/nmie-impl.cc
+        ${CMAKE_CURRENT_LIST_DIR}/pb11_wrapper.cc)
+
+# Define python extension
+add_library(python3-scattnlay SHARED ${_scattnlay_python_sources})
+target_link_libraries(python3-scattnlay PRIVATE Boost::headers ${Python_LIBRARIES})
+target_include_directories(python3-scattnlay PRIVATE ${NUMPY_INCLUDE_DIR} ${PYBIND11_INCLUDE_DIR})
+set_target_properties(python3-scattnlay PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
+
+set_target_properties(
+        python3-scattnlay
+        PROPERTIES
+        PREFIX ""
+        OUTPUT_NAME "scattnlay"
+        SUFFIX "${EXT_SUFFIX}"
+        LINKER_LANGUAGE C
+)
+
+# Sources for far field calculation
+set(_scattnlay_farfield_sources
+        ${CMAKE_CURRENT_LIST_DIR}/farfield.cc
+        ${CMAKE_CURRENT_LIST_DIR}/nmie.hpp
+        ${CMAKE_CURRENT_LIST_DIR}/nmie.cc)
+
+# Define exe for far field calculation
+add_executable(farfield ${_scattnlay_farfield_sources})
+target_link_libraries(farfield PRIVATE Boost::headers)
+set_target_properties(farfield PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
+
+# Sources for near field calculation
+set(_scattnlay_nearfield_sources
+        ${CMAKE_CURRENT_LIST_DIR}/nearfield.cc
+        ${CMAKE_CURRENT_LIST_DIR}/nmie.hpp
+        ${CMAKE_CURRENT_LIST_DIR}/nmie.cc)
+
+# Define exe for near field calculation
+add_executable(nearfield ${_scattnlay_nearfield_sources})
+target_link_libraries(nearfield PRIVATE Boost::headers)
+set_target_properties(nearfield PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
+
+# Rename files to match precision
+if (${ENABLE_MP})
+    set_property(TARGET python3-scattnlay APPEND_STRING PROPERTY OUTPUT_NAME "_mp")
+    set_property(TARGET farfield APPEND_STRING PROPERTY OUTPUT_NAME "farfield-mp")
+    set_property(TARGET nearfield APPEND_STRING PROPERTY OUTPUT_NAME "nearfield-mp")
+else ()
+    set_property(TARGET python3-scattnlay APPEND_STRING PROPERTY OUTPUT_NAME "_dp")
+    set_property(TARGET farfield APPEND_STRING PROPERTY OUTPUT_NAME "farfield-dp")
+    set_property(TARGET nearfield APPEND_STRING PROPERTY OUTPUT_NAME "nearfield-dp")
+endif ()

src/farfield.cc

@@ -29,17 +29,13 @@
 //    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
 //**********************************************************************************//
 
-#include <algorithm>
 #include <complex>
-#include <functional>
 #include <iostream>
 #include <stdexcept>
 #include <string>
 #include <vector>
-#include <stdlib.h>
-#include <stdio.h>
-#include <time.h>
-#include <string.h>
+#include <cstdio>
+
 #include "nmie.hpp"
 
 const double PI=3.14159265358979323846;
@@ -84,7 +80,7 @@ int main(int argc, char *argv[]) {
 
     int mode = -1;
     double tmp_mr;
-    for (auto arg : args) {
+    for (const auto& arg : args) {
       // For each arg in args list we detect the change of the current
       // read mode or read the arg. The reading args algorithm works
       // as a finite-state machine.
@@ -130,7 +126,7 @@ int main(int argc, char *argv[]) {
       }
 
       if (mode == read_mi) {
-        m.push_back(std::complex<double>( tmp_mr,std::stod(arg) ));
+        m.emplace_back( tmp_mr,std::stod(arg) );
         mode = read_x;
         continue;
       }
@@ -164,7 +160,7 @@ int main(int argc, char *argv[]) {
 
     if ( (x.size() != m.size()) || (L != x.size()) )
       throw std::invalid_argument(std::string("Broken structure!\n") + error_msg);
-    if ( (0 == m.size()) || ( 0 == x.size()) )
+    if ( (m.empty()) || ( x.empty()) )
       throw std::invalid_argument(std::string("Empty structure!\n") + error_msg);
 
     if (nt < 0) {

src/nearfield.cc

@@ -29,17 +29,13 @@
 //    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
 //**********************************************************************************//
 
-#include <algorithm>
 #include <complex>
-#include <functional>
 #include <iostream>
 #include <stdexcept>
 #include <string>
 #include <vector>
-#include <stdlib.h>
-#include <stdio.h>
-#include <time.h>
-#include <string.h>
+#include <cstdio>
+
 #include "nmie.hpp"
 
 const double PI=3.14159265358979323846;
@@ -84,7 +80,7 @@ int main(int argc, char *argv[]) {
 
     int mode = -1;
     double tmp_mr;
-    for (auto arg : args) {
+    for (const auto& arg : args) {
       // For each arg in args list we detect the change of the current
       // read mode or read the arg. The reading args algorithm works
       // as a finite-state machine.
@@ -130,7 +126,7 @@ int main(int argc, char *argv[]) {
       }
 
       if (mode == read_mi) {
-        m.push_back(std::complex<double>( tmp_mr,std::stod(arg) ));
+        m.emplace_back( tmp_mr,std::stod(arg) );
         mode = read_x;
         continue;
       }
@@ -211,7 +207,7 @@ int main(int argc, char *argv[]) {
 
     if ( (x.size() != m.size()) || (L != x.size()) )
       throw std::invalid_argument(std::string("Broken structure!\n") + error_msg);
-    if ( (0 == m.size()) || ( 0 == x.size()) )
+    if ( (m.empty()) || ( x.empty()) )
       throw std::invalid_argument(std::string("Empty structure!\n") + error_msg);
 
     if (nx == 1)

src/nmie-applied.cc

@@ -31,15 +31,12 @@
 //    @brief  Wrapper class around nMie function for ease of use                    //
 //                                                                                  //
 //**********************************************************************************//
+#include <stdexcept>
+#include <vector>
+
 #include "nmie-applied.hpp"
 #include "nmie-applied-impl.hpp"
 #include "nmie-precision.hpp"
-#include <array>
-#include <algorithm>
-#include <cstdio>
-#include <cstdlib>
-#include <stdexcept>
-#include <vector>
 
 namespace nmie {  
   //**********************************************************************************//
@@ -102,9 +99,7 @@ namespace nmie {
       // Will catch if  ml_mie fails or other errors.
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
-      return -1;
     }
-    return 0;
   }
   int nMieApplied(const unsigned int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
     return nmie::nMieApplied(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);

src/nmie-impl.cc

@@ -47,18 +47,22 @@
 //                                                                                  //
 // Hereinafter all equations numbers refer to [2]                                   //
 //**********************************************************************************//
-#include "nmie.hpp"
-#include "nmie-precision.hpp"
-#include <array>
-#include <algorithm>
-#include <cstdio>
-#include <cstdlib>
-#include <cmath>
+//<<<<<<< HEAD TODO
+//#include <array>
+//#include <algorithm>
+//#include <cstdio>
+//#include <cstdlib>
+//#include <cmath>
+//=======
+//>>>>>>> master
 #include <iostream>
 #include <iomanip>
 #include <stdexcept>
 #include <vector>
 
+#include "nmie.hpp"
+#include "nmie-precision.hpp"
+
 namespace nmie {
   //helper functions
 
@@ -1056,7 +1060,7 @@ namespace nmie {
 
     for (int n = nmax_ - 2; n >= 0; n--) {
       int n1 = n + 1;
-      FloatType rn = static_cast<FloatType>(n1);
+      auto rn = static_cast<FloatType>(n1);
 
       // using BH 4.12 and 4.50
       calcSpherHarm(Rho*ml, Theta, Phi, Psi[n1], D1n[n1], Pi[n], Tau[n], rn, M1o1n, M1e1n, N1o1n, N1e1n);

src/nmie-pybind11.cc

@@ -28,18 +28,11 @@
 //    You should have received a copy of the GNU General Public License             //
 //    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
 //**********************************************************************************//
-#include "nmie.hpp"
-#include "nmie-precision.hpp"
-#include <array>
-#include <tuple>
-#include <algorithm>
 #include <cstdio>
-#include <cstdlib>
-#include <stdexcept>
-#include <iostream>
-#include <iomanip>
 #include <vector>
 
+#include "nmie.hpp"
+
 #include <pybind11/pybind11.h>
 #include <pybind11/numpy.h>
 #include <pybind11/complex.h>
@@ -51,7 +44,7 @@ namespace py = pybind11;
 py::array_t< std::complex<double>> VectorComplex2Py(const std::vector<std::complex<double> > &c_x) {
   auto py_x = py::array_t< std::complex<double>>(c_x.size());
   auto py_x_buffer = py_x.request();
-  std::complex<double> *py_x_ptr = (std::complex<double> *) py_x_buffer.ptr;
+  auto *py_x_ptr = (std::complex<double> *) py_x_buffer.ptr;
   std::memcpy(py_x_ptr, c_x.data(), c_x.size()*sizeof(std::complex<double>));
   return py_x;
 }
@@ -72,11 +65,11 @@ std::vector<T> flatten(const std::vector<std::vector<T>>& v) {
 
 
 template <typename T>
-py::array VectorVector2Py(const std::vector<std::vector<T > > &x) {
+py::array Vector2DComplex2Py(const std::vector<std::vector<T > > &x) {
   size_t dim1 = x.size();
   size_t dim2 = x[0].size();
   auto result = flatten(x);
-  // https://github.com/tdegeus/pybind11_examples/blob/master/04_numpy-2D_cpp-vector/example.cpp 
+  // https://github.com/tdegeus/pybind11_examples/blob/master/04_numpy-2D_cpp-vector/example.cpp
   size_t              ndim    = 2;
   std::vector<size_t> shape   = {dim1, dim2};
   std::vector<size_t> strides = {sizeof(T)*dim2, sizeof(T)};
@@ -112,11 +105,29 @@ py::tuple py_ScattCoeffs(const py::array_t<double, py::array::c_style | py::arra
   std::vector<std::complex<double> > c_an, c_bn;
   int L = py_x.size();
   terms = nmie::ScattCoeffs(L, pl, c_x, c_m, nmax, c_an, c_bn);
-  
+
   return py::make_tuple(terms, VectorComplex2Py(c_an), VectorComplex2Py(c_bn));
 }
 
 
+py::tuple py_ExpanCoeffs(const py::array_t<double, py::array::c_style | py::array::forcecast> &py_x,
+                         const py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> &py_m,
+                         const int nmax=-1, const int pl=-1) {
+
+  auto c_x = Py2Vector<double>(py_x);
+  auto c_m = Py2Vector< std::complex<double> >(py_m);
+
+  int terms = 0;
+  int L = py_x.size();
+
+  std::vector<std::vector<std::complex<double> > > c_an(L + 1), c_bn(L + 1), c_cn(L + 1), c_dn(L + 1);
+
+  terms = nmie::ExpanCoeffs(L, pl, c_x, c_m, nmax, c_an, c_bn, c_cn, c_dn);
+
+  return py::make_tuple(terms, Vector2DComplex2Py(c_an), Vector2DComplex2Py(c_bn), Vector2DComplex2Py(c_cn), Vector2DComplex2Py(c_dn));
+}
+
+
 py::tuple py_scattnlay(const py::array_t<double, py::array::c_style | py::array::forcecast> &py_x,
                        const py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> &py_m,
                        const py::array_t<double, py::array::c_style | py::array::forcecast> &py_theta,
@@ -131,7 +142,7 @@ py::tuple py_scattnlay(const py::array_t<double, py::array::c_style | py::array:
   std::vector<std::complex<double> > c_S1, c_S2;
 
   terms = nmie::nMie(L, pl, c_x, c_m, nTheta, c_theta, nmax, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, c_S1, c_S2);
-  
+
   return py::make_tuple(terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo,
                         VectorComplex2Py(c_S1), VectorComplex2Py(c_S2));
 }
@@ -156,8 +167,8 @@ py::tuple py_fieldnlay(const py::array_t<double, py::array::c_style | py::array:
   for (auto& f : H) f.resize(3);
   int L = py_x.size(), terms;
   terms = nmie::nField(L, pl, c_x, c_m, nmax, nmie::Modes::kAll, nmie::Modes::kAll, ncoord, c_Xp, c_Yp, c_Zp, E, H);
-  auto py_E = VectorVector2Py<std::complex<double> >(E);
-  auto py_H = VectorVector2Py<std::complex<double> >(H);
+  auto py_E = Vector2DComplex2Py<std::complex<double> >(E);
+  auto py_H = Vector2DComplex2Py<std::complex<double> >(H);
   return py::make_tuple(terms, py_E, py_H);
 }
 

src/nmie-pybind11.hpp

@@ -43,6 +43,11 @@ py::tuple py_ScattCoeffs(const py::array_t<double, py::array::c_style | py::arra
                          const int nmax=-1, const int pl=-1);
 
 
+py::tuple py_ExpanCoeffs(const py::array_t<double, py::array::c_style | py::array::forcecast> &py_x,
+                         const py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> &py_m,
+                         const int nmax=-1, const int pl=-1);
+
+
 py::tuple py_scattnlay(const py::array_t<double, py::array::c_style | py::array::forcecast> &py_x,
                        const py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> &py_m,
                        const py::array_t<double, py::array::c_style | py::array::forcecast> &py_theta,

src/nmie.cc

@@ -45,20 +45,16 @@
 //                                                                                  //
 // Hereinafter all equations numbers refer to [2]                                   //
 //**********************************************************************************//
-#include "nmie.hpp"
-#include "nmie-precision.hpp"
-#include "nmie-impl.cc"
-#include <array>
-#include <algorithm>
-#include <cstdio>
-#include <cstdlib>
 #include <stdexcept>
 #include <iostream>
-#include <iomanip>
 #include <vector>
 
+#include "nmie.hpp"
+#include "nmie-precision.hpp"
+#include "nmie-impl.cc"
+
 namespace nmie {
-  
+
   //**********************************************************************************//
   // This function emulates a C call to calculate the scattering coefficients         //
   // required to calculate both the near- and far-field parameters.                   //
@@ -100,13 +96,11 @@ namespace nmie {
       // Will catch if  ml_mie fails or other errors.
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
-      return -1;
     }
-    return 0;
   }
 
   //**********************************************************************************//
-  // This function emulates a C call to calculate the scattering coefficients         //
+  // This function emulates a C call to calculate the expansion coefficients          //
   // required to calculate both the near- and far-field parameters.                   //
   //                                                                                  //
   // Input parameters:                                                                //
@@ -119,12 +113,14 @@ namespace nmie {
   //         set this parameter to -1 and the function will calculate it.             //
   //                                                                                  //
   // Output parameters:                                                               //
-  //   aln, bln ,cln, dln : Complex scattering amplitudes                             //
+  //   aln, bln, cln, dln: Complex expansion coefficients                             //
   //                                                                                  //
   // Return value:                                                                    //
-  //   Number of multipolar expansion terms used for the calculations in each layer   //
+  //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int ExpansionCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, std::vector<std::vector<std::complex<double> > >& aln, std::vector<std::vector<std::complex<double> > >& bln, std::vector<std::vector<std::complex<double> > >& cln, std::vector<std::vector<std::complex<double> > >& dln) {
+  int ExpanCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax,
+                  std::vector<std::vector<std::complex<double> > >& aln, std::vector<std::vector<std::complex<double> > >& bln,
+                  std::vector<std::vector<std::complex<double> > >& cln, std::vector<std::vector<std::complex<double> > >& dln) {
 
     if (x.size() != L || m.size() != L)
         throw std::invalid_argument("Declared number of layers do not fit x and m!");
@@ -135,24 +131,25 @@ namespace nmie {
       ml_mie.SetPECLayer(pl);
       ml_mie.SetMaxTerms(nmax);
 
+    // Calculate scattering coefficients an_ and bn_
       ml_mie.calcScattCoeffs();
+    // Calculate expansion coefficients aln_,  bln_, cln_, and dln_
       ml_mie.calcExpanCoeffs();
 
-      aln = ConvertComplexVectorVector<double>(ml_mie.GetAln());
-      bln = ConvertComplexVectorVector<double>(ml_mie.GetBln());
-      cln = ConvertComplexVectorVector<double>(ml_mie.GetCln());
-      dln = ConvertComplexVectorVector<double>(ml_mie.GetDln());
+      aln = ConvertComplexVectorVector<double>(ml_mie.GetLayerAn());
+      bln = ConvertComplexVectorVector<double>(ml_mie.GetLayerBn());
+      cln = ConvertComplexVectorVector<double>(ml_mie.GetLayerCn());
+      dln = ConvertComplexVectorVector<double>(ml_mie.GetLayerDn());
 
       return ml_mie.GetMaxTerms();
     } catch(const std::invalid_argument& ia) {
       // Will catch if  ml_mie fails or other errors.
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
-      return -1;
     }
-    return 0;
   }
 
+
   //**********************************************************************************//
   // This function emulates a C call to calculate the actual scattering parameters    //
   // and amplitudes.                                                                  //
@@ -208,7 +205,7 @@ namespace nmie {
       // 	<< "Qext = "
       // 	<< ml_mie.GetQext()
       // 	<< std::endl;
-      
+
       *Qext = static_cast<double>(ml_mie.GetQext());
       *Qsca = static_cast<double>(ml_mie.GetQsca());
       *Qabs = static_cast<double>(ml_mie.GetQabs());
@@ -224,9 +221,7 @@ namespace nmie {
       // Will catch if  ml_mie fails or other errors.
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
-      return -1;
     }
-    return 0;
   }
 
 
@@ -404,10 +399,10 @@ namespace nmie {
     if (Xp_vec.size() != ncoord || Yp_vec.size() != ncoord || Zp_vec.size() != ncoord
         || E.size() != ncoord || H.size() != ncoord)
       throw std::invalid_argument("Declared number of coords do not fit Xp, Yp, Zp, E, or H!");
-    for (auto f:E)
+    for (const auto& f:E)
       if (f.size() != 3)
         throw std::invalid_argument("Field E is not 3D!");
-    for (auto f:H)
+    for (const auto& f:H)
       if (f.size() != 3)
         throw std::invalid_argument("Field H is not 3D!");
     try {
@@ -430,8 +425,6 @@ namespace nmie {
       // Will catch if  ml_mie fails or other errors.
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
-      return - 1;
     }
-    return 0;
   }
 }  // end of namespace nmie

src/nmie.hpp

@@ -42,9 +42,13 @@
 #endif
 namespace nmie {
   int ScattCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn);
-  int ExpansionCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, std::vector<std::vector<std::complex<double> > >& an, std::vector<std::vector<std::complex<double> > >& bn, std::vector<std::vector<std::complex<double> > >& cn, std::vector<std::vector<std::complex<double> > >& dn);
+
+  int ExpanCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax,
+                std::vector<std::vector<std::complex<double> > >& an, std::vector<std::vector<std::complex<double> > >& bn,
+                std::vector<std::vector<std::complex<double> > >& cn, std::vector<std::vector<std::complex<double> > >& dn) {
+
   // pl, nmax, mode_n, mode_type
-  int nMie(const unsigned int L,
+    int nMie(const unsigned int L,
            const int pl,
            std::vector<double>& x, std::vector<std::complex<double> >& m,
            const unsigned int nTheta, std::vector<double>& Theta,
@@ -54,7 +58,7 @@ namespace nmie {
            std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2,
            int mode_n, int mode_type);
   // pl and nmax
-  int nMie(const unsigned int L,
+    int nMie(const unsigned int L,
            const int pl,
            std::vector<double>& x, std::vector<std::complex<double> >& m,
            const unsigned int nTheta, std::vector<double>& Theta,
@@ -85,18 +89,17 @@ namespace nmie {
            double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr,
            double *g, double *Albedo,
            std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
-  int nField(const unsigned int L, const int pl, const std::vector<double>& x,
-          const std::vector<std::complex<double> >& m, const int nmax,
-          const int mode_n, const int mode_type,
-          const unsigned int ncoord,
-          const std::vector<double>& Xp, const std::vector<double>& Yp, const std::vector<double>& Zp,
-          std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H);
+  int nField(const unsigned int L, const int pl,
+             const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax,
+             const int mode_n, const int mode_type, const unsigned int ncoord,
+             const std::vector<double>& Xp, const std::vector<double>& Yp, const std::vector<double>& Zp,
+             std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H);
 
   // constants for per mode evaluation
   enum Modes {kAll = -1, kElectric = 0, kMagnetic = 1};
 
   template <typename FloatType = double>
-  class MultiLayerMie {    
+  class MultiLayerMie {
    public:
     //Used constants TODO! Change to boost PI
     const double PI_=3.14159265358979323846;
@@ -124,10 +127,10 @@ namespace nmie {
     std::vector<std::complex<FloatType> > GetAn(){return an_;};
     std::vector<std::complex<FloatType> > GetBn(){return bn_;};
 
-    std::vector<std::vector<std::complex<FloatType> > > GetAln(){return aln_;};
-    std::vector<std::vector<std::complex<FloatType> > > GetBln(){return bln_;};
-    std::vector<std::vector<std::complex<FloatType> > > GetCln(){return cln_;};
-    std::vector<std::vector<std::complex<FloatType> > > GetDln(){return dln_;};
+    std::vector< std::vector<std::complex<FloatType> > > GetLayerAn(){return aln_;};
+    std::vector< std::vector<std::complex<FloatType> > > GetLayerBn(){return bln_;};
+    std::vector< std::vector<std::complex<FloatType> > > GetLayerCn(){return cln_;};
+    std::vector< std::vector<std::complex<FloatType> > > GetLayerDn(){return dln_;};
 
     // Problem definition
     // Modify size of all layers
@@ -209,7 +212,7 @@ namespace nmie {
     void calcSpherHarm(const std::complex<FloatType> Rho, const FloatType Theta, const FloatType Phi,
                        const std::complex<FloatType>& rn, const std::complex<FloatType>& Dn,
                        const FloatType& Pi, const FloatType& Tau, const FloatType& n,
-                       std::vector<std::complex<FloatType> >& Mo1n, std::vector<std::complex<FloatType> >& Me1n, 
+                       std::vector<std::complex<FloatType> >& Mo1n, std::vector<std::complex<FloatType> >& Me1n,
                        std::vector<std::complex<FloatType> >& No1n, std::vector<std::complex<FloatType> >& Ne1n);
 
     void calcFieldByComponents(const FloatType Rho, const FloatType Theta, const FloatType Phi,

src/pb11_wrapper.cc

@@ -16,6 +16,10 @@ PYBIND11_MODULE(scattnlay_dp, m)
         "Calculate the scattering coefficients, required to calculate both the near- and far-field parameters.",
         py::arg("x"), py::arg("m"), py::arg("nmax") = -1, py::arg("pl") = -1);
 
+  m.def("expancoeffs", &py_ExpanCoeffs,
+        "Calculate the expansion coefficients, required to calculate the near-field parameters.",
+        py::arg("x"), py::arg("m"), py::arg("nmax") = -1, py::arg("pl") = -1);
+
   m.def("scattnlay", &py_scattnlay,
         "Calculate the scattering parameters and amplitudes.",
         py::arg("x"), py::arg("m"), py::arg("theta") = py::array_t<double>(0), py::arg("nmax") = -1, py::arg("pl") = -1);

src/shell-generator.cc

@@ -30,18 +30,15 @@
 //**********************************************************************************//
 //   @brief  Generates points for integration on sphere surface
 
-#include <algorithm>
 #include <cmath>
 #include <complex>
-#include <cstdio>
-#include <fstream>
 #include <iostream>
 #include <set>
-#include <sstream>
 #include <stdexcept>
-#include <string>
 #include <vector>
+
 #include "shell-generator.hpp"
+
 namespace shell_generator {
   template<class T> inline T pow2(const T value) {return value*value;}
   // ********************************************************************** //