new poynting plotting

tests/python/calc-SiAgSi-Qabs.py

@@ -0,0 +1,277 @@
+#!/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.
+
+import scattnlay
+from scattnlay import fieldnlay
+from scattnlay import scattnlay
+import numpy as np
+import cmath
+from fieldplot import GetFlow3D
+from fieldplot import GetField
+
+###############################################################################
+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
+design = 5   # Bulk Ag
+x, m, WL = SetXM(design)
+
+
+WL_units='nm'
+comment='P-SiAgSi-flow'
+comment='bulk-P-Ag-flow'
+print "x =", x
+print "m =", m
+npts = 101
+factor=2.2
+flow_total = 3
+#flow_total = 0
+crossplane='XZ'
+#crossplane='YZ'
+#crossplane='XY'
+
+
+
+Ec, Hc, P, coordX, coordZ = GetField(crossplane, npts, factor, x, m)
+
+Er = np.absolute(Ec)
+Hr = np.absolute(Hc)
+# |E|/|Eo|
+Eabs = np.sqrt(Er[ :, 0]**2 + Er[ :, 1]**2 + Er[ :, 2]**2)
+Eangle = np.angle(Ec[ :, 0])/np.pi*180
+Habs= np.sqrt(Hr[ :, 0]**2 + Hr[ :, 1]**2 + Hr[ :, 2]**2)
+Hangle = np.angle(Hc[ :, 1])/np.pi*180
+
+try:
+    import matplotlib.pyplot as plt
+    from matplotlib import cm
+    from matplotlib.colors import LogNorm
+
+    Eabs_data = np.resize(P, (npts, npts)).T
+    #Eabs_data = np.resize(Pabs, (npts, npts)).T
+    # Eangle_data = np.resize(Eangle, (npts, npts)).T
+    # Habs_data = np.resize(Habs, (npts, npts)).T
+    # Hangle_data = np.resize(Hangle, (npts, npts)).T
+
+    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('Pabs')
+    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=1, 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]*15
+        #max_length=factor*x[-1]*5
+        max_angle = np.pi/200
+        #for flow in range(0,flow_total*2+1):
+        for flow in range(0,flow_total):
+            if crossplane=='XZ':
+                #x0 = min_SP*2 + flow*step_SP
+                x0 = min_SP + flow*step_SP
+                z0 = min_SP
+                #y0 = x[-1]/20 
+            elif crossplane=='YZ':
+                #y0 = min_SP*2 + flow*step_SP
+                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)
+            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=1, 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+".svg")
+    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/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
 #
@@ -137,7 +138,7 @@ coord = np.vstack((coordX, coordY, coordZ)).transpose()
 
 terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
 terms, E, H = fieldnlay(x, m, coord)
-
+print(E)
 P = np.array(map(lambda n: np.linalg.norm(np.cross(E[0][n], H[0][n])).real, range(0, len(E[0]))))
 
 Ec = E[0, :, :]

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,168 @@
+#!/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/>.
+
+# 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):
+    # 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)
+    Ec, Hc = E[0, 0, :], H[0, 0, :]
+    S = np.cross(Ec, Hc.conjugate()).real
+    Snorm_prev = S/np.linalg.norm(S)
+    length = 0
+    dpos = step
+    while length < max_length:
+        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 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)
+            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
+            Snorm = S/np.linalg.norm(S)
+            #Snorm = Snorm.real
+            #Snorm = np.absolute(Snorm)
+            diff = Snorm-Snorm_prev
+            if np.linalg.norm(diff)<0.05:
+            # angle = angle_between(Snorm, Snorm_prev)
+            # if abs(angle) < max_angle:
+                break
+            step = step/2.0
+        #3. Save result
+        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):
+    """
+    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)
+    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
+
+