Merge branch 'master' of github.com:ovidiopr/scattnlay

Conflicts:
	tests/python/field-Ag-flow.py

nmie.cc

@@ -1144,7 +1144,7 @@ namespace nmie {
       if (std::abs(bln_[0][n]) < 1e-10) bln_[0][n] = 0.0;
       else {
         //throw std::invalid_argument("Unstable calculation of bln_[0][n]!");
-        printf("Warning: Potentially unstable calculation of bln (bln[0][%i] = %g, %gi)\n", n, bln_[0][n].real(), bln_[0][n].imag());
+        printf("Warning: Potentially unstable calculation of bln (bln[0][%i] = %g, %gi) pl=%d\n", n, bln_[0][n].real(), bln_[0][n].imag(), PEC_layer_position_);
         bln_[0][n] = 0.0;
       }
     }

tests/python/calc-SiAgSi-Qabs.py

@@ -0,0 +1,153 @@
+#!/usr/bin/env python
+# -*- coding: UTF-8 -*-
+#
+#    Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
+#    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>
+#
+#    This file is part of python-scattnlay
+#
+#    This program is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation, either version 3 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    The only additional remark is that we expect that all publications
+#    describing work using this software, or all commercial products
+#    using it, cite the following reference:
+#    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
+#        a multilayered sphere," Computer Physics Communications,
+#        vol. 180, Nov. 2009, pp. 2348-2354.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+# This test case calculates the electric field in the 
+# E-k plane, for an spherical Si-Ag-Si nanoparticle.
+
+import scattnlay
+from scattnlay import fieldnlay
+from scattnlay import scattnlay
+import numpy as np
+import cmath
+# from fieldplot import GetFlow3D
+# from fieldplot import GetField
+from fieldplot import fieldplot
+
+###############################################################################
+def SetXM(design):
+    """ design value:
+    1: AgSi - a1
+    2: SiAgSi - a1, b1
+    3: SiAgSi - a1, b2
+    """
+    epsilon_Si = 18.4631066585 + 0.6259727805j
+    epsilon_Ag = -8.5014154589 + 0.7585845411j
+    index_Si = np.sqrt(epsilon_Si)
+    index_Ag = np.sqrt(epsilon_Ag)
+    isSiAgSi=True
+    isBulk = False
+    if design==1:
+        #36	5.62055	0	31.93	4.06	49	5.62055	500
+        isSiAgSi=False
+        WL=500 #nm
+        core_width = 0.0 #nm Si
+        inner_width = 31.93 #nm Ag
+        outer_width = 4.06 #nm  Si
+    elif design==2:
+        #62.5	4.48866	29.44	10.33	22.73	0	4.48866	500
+        WL=500 #nm
+        core_width = 29.44 #nm Si
+        inner_width = 10.33 #nm Ag
+        outer_width = 22.73 #nm  Si
+    elif design == 3:
+        #81.4	3.14156	5.27	8.22	67.91	0	3.14156	500
+        WL=500 #nm
+        core_width = 5.27 #nm Si
+        inner_width = 8.22 #nm Ag
+        outer_width = 67.91 #nm  Si
+    elif design==4:
+        WL=800 #nm
+        epsilon_Si = 13.64 + 0.047j
+        epsilon_Ag = -28.05 + 1.525j
+        core_width = 17.74 #nm Si
+        inner_width = 23.31 #nm Ag
+        outer_width = 22.95 #nm  Si
+    elif design==5:
+        WL=354 #nm
+        core_r = WL/20.0
+        epsilon_Ag = -2.0 + 0.28j   #original
+        index_Ag = np.sqrt(epsilon_Ag)
+        x = np.ones((1), dtype = np.float64)
+        x[0] = 2.0*np.pi*core_r/WL
+        m = np.ones((1), dtype = np.complex128)
+        m[0] = index_Ag
+        # x = np.ones((2), dtype = np.float64)
+        # x[0] = 2.0*np.pi*core_r/WL/4.0*3.0
+        # x[1] = 2.0*np.pi*core_r/WL
+        # m = np.ones((2), dtype = np.complex128)
+        # m[0] = index_Ag
+        # m[1] = index_Ag
+        return x, m, WL
+
+
+    core_r = core_width
+    inner_r = core_r+inner_width
+    outer_r = inner_r+outer_width
+
+    nm = 1.0
+    if isSiAgSi:
+        x = np.ones((3), dtype = np.float64)
+        x[0] = 2.0*np.pi*core_r/WL
+        x[1] = 2.0*np.pi*inner_r/WL
+        x[2] = 2.0*np.pi*outer_r/WL
+        m = np.ones((3), dtype = np.complex128)
+        m[0] = index_Si/nm
+        m[1] = index_Ag/nm
+    #    m[0, 1] = index_Si/nm
+        m[2] = index_Si/nm
+    else:
+        # bilayer
+        x = np.ones((2), dtype = np.float64)
+        x[0] = 2.0*np.pi*inner_r/WL
+        x[1] = 2.0*np.pi*outer_r/WL
+        m = np.ones((2), dtype = np.complex128)
+        m[0] = index_Ag/nm
+        m[1] = index_Si/nm
+    return x, m, WL
+
+
+###############################################################################
+#design = 1 #AgSi
+design = 2
+#design = 3
+#design = 4   # WL=800
+comment='SiAgSi-flow'
+#design = 5   # Bulk Ag
+# comment='bulk-Ag-flow'
+x, m, WL = SetXM(design)
+
+WL_units='nm'
+print "x =", x
+print "m =", m
+npts = 501
+factor=2.1
+flow_total = 39
+#flow_total = 21
+#flow_total = 0
+crossplane='XZ'
+#crossplane='YZ'
+#crossplane='XY'
+
+# Options to plot: Eabs, Habs, Pabs, angleEx, angleHy
+field_to_plot='Pabs'
+#field_to_plot='angleEx'
+
+fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total)
+
+
+

tests/python/coating-flow.py

@@ -0,0 +1,121 @@
+#!/usr/bin/env python
+# -*- coding: UTF-8 -*-
+#
+#    Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
+#
+#    This file is part of python-scattnlay
+#
+#    This program is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation, either version 3 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    The only additional remark is that we expect that all publications
+#    describing work using this software, or all commercial products
+#    using it, cite the following reference:
+#    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
+#        a multilayered sphere," Computer Physics Communications,
+#        vol. 180, Nov. 2009, pp. 2348-2354.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+# This test case calculates the electric field in the
+# E-k plane, for an spherical Si-Ag-Si nanoparticle. Core radius is 17.74 nm,
+# inner layer 23.31nm, outer layer 22.95nm. Working wavelength is 800nm, we use
+# silicon epsilon=13.64+i0.047, silver epsilon= -28.05+i1.525
+
+import os, cmath
+import numpy as np
+from scattnlay import fieldnlay
+from fieldplot import fieldplot
+
+
+if __name__ == '__main__':
+    import argparse
+
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument("dirnames", nargs='*', default='.', help="read all data from DIR(S)")
+    parser.add_argument("-f", "--filename", dest="fname", nargs='?', default=None,
+                        help="name of 'n' file")
+    parser.add_argument("-w", "--wavelength", dest="wl", default=3.75, type=float,
+                        help="wavelength of electromagnetic wave")
+    parser.add_argument("-r", "--radius", dest="rad", default=None, type=float,
+                        help="radius of PEC sphere")
+    parser.add_argument("-t", "--thickness", dest="tc", default=0.8, type=float,
+                        help="thickness of cloaking layer")
+    parser.add_argument("-n", "--npoints", dest="npts", default=101, type=int,
+                        help="number of points for the grid")
+
+    args = parser.parse_args()
+
+    for dirname in args.dirnames:
+        print "Calculating spectra for data file(s) in dir '%s'..." % (dirname)
+
+        wl = args.wl # cm
+        if (args.rad is None):
+            Rs = 0.75*wl  # cm
+        else:
+            Rs = args.rad # cm
+        tc = args.tc # cm
+
+        if (args.fname is None):
+            files = [x for x in os.listdir('%s/' % (dirname)) if x.endswith('.dat')]
+            files.sort()
+        else:
+            files = [args.fname]
+
+        npts = args.npts # cm
+
+        if not os.path.exists('%s/flow-results/' % (dirname)):
+            os.makedirs('%s/flow-results/' % (dirname))
+
+        Rt = Rs + tc # cm
+
+        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)
+
+            basename = os.path.splitext(fname)[0]
+
+            nvalues = np.loadtxt('%s/%s' % (dirname, fname))*1.0 + 1e-11j
+
+            tl = tc/len(nvalues) # cm
+            r = [Rs]
+
+            for i in range(len(nvalues)):
+                r += [r[i] + tl]
+
+            x = np.ones((1, len(nvalues) + 1), dtype = np.float64)
+            m = np.ones((1, len(nvalues) + 1), dtype = np.complex128)
+
+            x[0] = 2.0*np.pi*np.array(r, dtype = np.float64)/wl
+            m[0] = np.array([ms] + nvalues[:, 1].tolist(), dtype = np.complex128)
+
+            factor = 2
+            comment='PEC-'+basename
+            WL_units='cm'
+            flow_total = 39
+            #flow_total = 0
+            #crossplane='XZ'
+            crossplane='YZ'
+            #crossplane='XY'
+
+            # Options to plot: Eabs, Habs, Pabs, angleEx, angleHy
+            field_to_plot='Pabs'
+            #field_to_plot='Eabs'
+            #field_to_plot='angleEx'
+            fieldplot(x[0],m[0], wl, comment, WL_units, crossplane, field_to_plot, npts,
+                      factor, flow_total, pl=0, outline_width=0.1)
+
+
+        print "Done!!"
+

tests/python/coating-flow.sh

@@ -0,0 +1,5 @@
+#!/bin/sh
+
+python coating-flow.py -w 3.75 -t 0.80 -f index-dv.dat -n 1501
+#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

tests/python/field-Ag-flow.py

@@ -2,6 +2,7 @@
 # -*- coding: UTF-8 -*-
 #
 #    Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
+#    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>
 #
 #    This file is part of python-scattnlay
 #
@@ -26,16 +27,17 @@
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 # This test case calculates the electric field in the 
-# E-k plane, for an spherical Si-Ag-Si nanoparticle. Core radius is 17.74 nm,
-# inner layer 23.31nm, outer layer 22.95nm. Working wavelength is 800nm, we use
-# silicon epsilon=13.64+i0.047, silver epsilon= -28.05+i1.525
+# E-k plane, for an spherical Ag nanoparticle.
 
 import scattnlay
 from scattnlay import fieldnlay
 from scattnlay import scattnlay
+from fieldplot import fieldplot
+
 import numpy as np
 import cmath
 
+<<<<<<< HEAD
 
 def get_index(array,value):
     idx = (np.abs(array-value)).argmin()
@@ -80,6 +82,8 @@ def GetFlow(scale_x, scale_z, Ec, Hc, a, b, nmax):
         flow_z.append(z_pos)
     return flow_x, flow_z
 
+=======
+>>>>>>> feb3ad9a4b3aa424f2e1087b4bc7b9bc52598810
 # # a)
 #WL=400 #nm
 #core_r = WL/20.0
@@ -100,25 +104,23 @@ epsilon_Ag = -2.0 + 0.28j
 #core_r = WL/20.0
 #epsilon_Ag = -2.71 + 0.25j
 
-
 index_Ag = np.sqrt(epsilon_Ag)
-
-
 # n1 = 1.53413
 # n2 = 0.565838 + 7.23262j
 nm = 1.0
 
-x = np.ones((1, 2), dtype = np.float64)
-x[0, 0] = 2.0*np.pi*core_r/WL/4.0*3.0
-x[0, 1] = 2.0*np.pi*core_r/WL
+x = np.ones((2), dtype = np.float64)
+x[0] = 2.0*np.pi*core_r/WL/4.0*3.0
+x[1] = 2.0*np.pi*core_r/WL
 
-m = np.ones((1, 2), dtype = np.complex128)
-m[0, 0] = index_Ag/nm
-m[0, 1] = index_Ag/nm
+m = np.ones((2), dtype = np.complex128)
+m[0] = index_Ag/nm
+m[1] = index_Ag/nm
 
 print "x =", x
 print "m =", m
 
+<<<<<<< HEAD
 npts = 281
 
 factor=3
@@ -248,5 +250,24 @@ try:
 finally:
     print("Qabs = "+str(Qabs));
 #
+=======
+comment='bulk-Ag-flow'
+WL_units='nm'
+npts = 501
+factor=2.1
+flow_total = 39
+#flow_total = 21
+#flow_total = 0
+crossplane='XZ'
+#crossplane='YZ'
+#crossplane='XY'
+
+# Options to plot: Eabs, Habs, Pabs, angleEx, angleHy
+field_to_plot='Pabs'
+#field_to_plot='angleEx'
+
+fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total, is_flow_extend=False)
+
+>>>>>>> feb3ad9a4b3aa424f2e1087b4bc7b9bc52598810
 
 

tests/python/field-SiAgSi.py

@@ -2,6 +2,7 @@
 # -*- coding: UTF-8 -*-
 #
 #    Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
+#    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>
 #
 #    This file is part of python-scattnlay
 #

tests/python/fieldplot.py

@@ -0,0 +1,320 @@
+#!/usr/bin/env python
+# -*- coding: UTF-8 -*-
+#
+#    Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
+#    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>
+#
+#    This file is part of python-scattnlay
+#
+#    This program is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation, either version 3 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    The only additional remark is that we expect that all publications
+#    describing work using this software, or all commercial products
+#    using it, cite the following reference:
+#    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
+#        a multilayered sphere," Computer Physics Communications,
+#        vol. 180, Nov. 2009, pp. 2348-2354.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+# Several functions to plot field and streamlines (power flow lines).
+
+import scattnlay
+from scattnlay import fieldnlay
+from scattnlay import 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
+    coord = np.vstack(([flow_x[-1]], [flow_y[-1]], [flow_z[-1]])).transpose()
+    terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord,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>3000): # 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(([flow_x[-1]+dx], [flow_y[-1]+dy], [flow_z[-1]+dz])).transpose()
+            terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord,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 xrange(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):
+    """
+    crossplane: XZ, YZ, XY
+    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)
+        coordY = zero
+        coordPlot1 = coordX
+        coordPlot2 = coordZ
+    elif crossplane=='YZ':
+        coordY, coordZ = np.meshgrid(scan, scan)
+        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)
+        coordZ = zero
+        coordPlot1 = coordY
+        coordPlot2 = coordX
+        
+    coord = np.vstack((coordX, coordY, coordZ)).transpose()
+    terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord, pl=pl)
+    Ec = E[0, :, :]
+    Hc = H[0, :, :]
+    P=[]
+    P = np.array(map(lambda n: np.linalg.norm(np.cross(Ec[n], Hc[n])).real, range(0, len(E[0]))))
+
+    # 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
+###############################################################################
+def fieldplot(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):
+    Ec, Hc, P, coordX, coordZ = GetField(crossplane, npts, factor, x, m, pl)
+    Er = np.absolute(Ec)
+    Hr = np.absolute(Hc)
+    try:
+        import matplotlib.pyplot as plt
+        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)$'
+        elif field_to_plot == 'Eabs':
+            Eabs = np.sqrt(Er[ :, 0]**2 + Er[ :, 1]**2 + Er[ :, 2]**2)
+            Eabs_data = np.resize(Eabs, (npts, npts)).T 
+            label = r'$|E|$'
+        elif field_to_plot == 'Habs':
+            Habs= np.sqrt(Hr[ :, 0]**2 + Hr[ :, 1]**2 + Hr[ :, 2]**2)
+            Eabs_data = np.resize(Habs, (npts, npts)).T 
+            label = r'$|H|$'
+        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)$'
+        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)$'
+
+        fig, ax = plt.subplots(1,1)
+        # 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)
+
+        # Define scale ticks
+        min_tick = np.amin(Eabs_data[~np.isnan(Eabs_data)])
+        max_tick = np.amax(Eabs_data[~np.isnan(Eabs_data)])
+        scale_ticks = np.linspace(min_tick, max_tick, 6)
+
+        # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
+        ax.set_title(label)
+        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))
+                        #,norm = LogNorm()
+                        )
+        ax.axis("image")
+
+        # Add colorbar
+        cbar = fig.colorbar(cax, ticks = [a for a in scale_ticks])
+        cbar.ax.set_yticklabels(['%5.3g' % (a) for a in scale_ticks]) # vertically oriented colorbar
+        pos = list(cbar.ax.get_position().bounds)
+        #fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
+        if crossplane=='XZ':
+            plt.xlabel('Z, '+WL_units)
+            plt.ylabel('X, '+WL_units)
+        elif crossplane=='YZ':
+            plt.xlabel('Z, '+WL_units)
+            plt.ylabel('Y, '+WL_units)
+        elif crossplane=='XY':
+            plt.xlabel('Y, '+WL_units)
+            plt.ylabel('X, '+WL_units)
+
+
+        # # 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 (crossplane=='XZ' or crossplane=='YZ') and flow_total>0:
+
+            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]*8
+            #max_length=factor*x[-1]*5
+            max_angle = np.pi/160
+            if is_flow_extend:
+                rg = range(0,flow_total*2+1)
+            else:
+                rg = range(0,flow_total)
+            for flow in rg:
+                if crossplane=='XZ':
+                    if is_flow_extend:
+                        x0 = min_SP*2 + flow*step_SP
+                    else:
+                        x0 = min_SP + flow*step_SP
+                    z0 = min_SP
+                    #y0 = x[-1]/20 
+                elif crossplane=='YZ':
+                    if is_flow_extend:
+                        y0 = min_SP*2 + flow*step_SP
+                    else:
+                        y0 = min_SP + flow*step_SP
+                    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
+
+                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=0.7, edgecolor='white',zorder = 1.9)
+                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')
+
+        plt.savefig(comment+"-R"+str(int(round(x[-1]*WL/2.0/np.pi)))+"-"+crossplane+"-"
+                    +field_to_plot+".pdf")
+        plt.draw()
+
+    #    plt.show()
+
+        plt.clf()
+        plt.close()
+    finally:
+        terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(np.array([x]),
+                                                                         np.array([m]))
+        print("Qabs = "+str(Qabs));
+    #
+
+

tests/python/index-ch.dat

@@ -0,0 +1,98 @@
+#Layer N	Index	
+# TargetR2.8125-index-n-1k-1-CoatingW02.400-FinalRCS24.0577diff-54.3%-n96-s406.642299335648-index
+1	5.749335803512620302	
+2	3.141762961744382743	
+3	3.517942484287340577	
+4	4.791493884037310025	
+5	3.735822254252520125	
+6	2.083879849088456027	
+7	6.690275153118312446	
+8	4.93218340173765224	
+9	1.775615487238558776	
+10	5.176026567191519767	
+11	3.743194574611140357	
+12	2.866197219283298292	
+13	2.185180285910699904	
+14	4.414273003084785429	
+15	4.868023736776916799	
+16	1.710170432391318052	
+17	3.50796206700414448	
+18	2.697988996898661895	
+19	3.553957623718730474	
+20	7.999980780105266476	
+21	3.894242708766762373	
+22	7.998130485037586723	
+23	2.331279488473009032	
+24	6.079118361707619655	
+25	7.432492661618770313	
+26	4.595789892550891764	
+27	4.292344676729104691	
+28	5.332981371688163108	
+29	2.93942304734845461	
+30	4.706237624636468553	
+31	3.348093523938082239	
+32	1.725161528985504011	
+33	2.700514574998721606	
+34	2.749975835488200016	
+35	6.127403895715513471	
+36	7.999989501592007279	
+37	7.999999862141754292	
+38	7.997222358772536843	
+39	7.885835231178068838	
+40	1.104874591509938497	
+41	2.557385992494067839	
+42	1.283487106436893566	
+43	4.280402167755322118	
+44	4.15545727580456159	
+45	1.107605812272585144	
+46	1.66302549303807945	
+47	1.406567313133828545	
+48	1.638577522082844684	
+49	1.40392641764507875	
+50	6.578845533119020317	
+51	3.58680531520020951	
+52	4.162243508120385904	
+53	2.420274012058465374	
+54	3.806898234367299239	
+55	7.9999154826863208	
+56	7.299283239081675312	
+57	7.957581641865791866	
+58	7.999756463939247553	
+59	1.000307210296912386	
+60	1.268650135088047692	
+61	6.558667066646631127	
+62	1.009636547370579063	
+63	7.683437169416794887	
+64	7.616087731862954513	
+65	1.412789224110302388	
+66	5.092962017268542141	
+67	2.636203543845038411	
+68	5.864195462731800923	
+69	5.559535453969105312	
+70	4.212327947425384522	
+71	3.775995746343467463	
+72	2.79644746238645947	
+73	2.990656024028564897	
+74	2.268934872808900938	
+75	3.53529066398966263	
+76	5.170244446423582829	
+77	1.000153502850401388	
+78	4.146432814237805964	
+79	4.502227612310949922	
+80	3.574165620116289244	
+81	5.941097162753926852	
+82	7.999414796963508323	
+83	5.936878282510322791	
+84	4.649745206189950686	
+85	7.572390563536975705	
+86	4.780985022403201334	
+87	4.356413030716980828	
+88	3.718350395661014218	
+89	4.11922953620705723	
+90	3.151317818732550791	
+91	2.602738005170567526	
+92	4.780286744896858231	
+93	1.297975845438861464	
+94	6.370021149298636942	
+95	1.418719383122659439	
+96	4.580971796543868102	

tests/python/index-dv.dat

@@ -0,0 +1,34 @@
+#Layer N	Index	
+# TargetR2.8125-CoatingW00.800-FinalRCS25.3412diff-51.9%-n32-s139.585414476121-index
+1	2.914401365946849953	
+2	7.999999568066106903	
+3	7.999999805091645655	
+4	7.999999238794624823	
+5	7.999999854193549531	
+6	5.979617734581032629	
+7	2.90056030279753907	
+8	2.513877134657918955	
+9	2.313975248489835934	
+10	2.186263097301196279	
+11	2.103731626279280675	
+12	2.050425337776427881	
+13	2.02722718503005872	
+14	2.043407875278105301	
+15	2.141954143805122612	
+16	2.779546399372911925	
+17	7.999999963405493908	
+18	7.999999959245399417	
+19	7.999999992785230774	
+20	7.999999763204644232	
+21	7.999999962594151803	
+22	3.607211241982554384	
+23	2.519963615602491824	
+24	2.30033851793211408	
+25	2.210019176974556299	
+26	2.177346881906367049	
+27	2.187732694221404373	
+28	2.241437964316892462	
+29	2.363784302705050688	
+30	2.652539546222675515	
+31	5.370055597183615248	
+32	7.999999378375875381	

tests/python/index-sv.dat

@@ -0,0 +1,18 @@
+#Layer N	Index	
+# TargetR2.8125-CoatingW00.620-FinalRCS27.4228diff-47.9%-n16-s50.421839621567-index
+1	1.000000000922160126	
+2	1.000000697299052721	
+3	7.994698322241974076	
+4	7.999985996968639768	
+5	6.748281789975593803	
+6	3.482762722162508151	
+7	2.380626616473425994	
+8	2.159309340071060568	
+9	2.034109569517063676	
+10	1.884925694693077025	
+11	1.828250073364771744	
+12	2.04069028461526214	
+13	2.296151028318206055	
+14	5.572006973799219054	
+15	1.000005309961293776	
+16	1.000035201183201217