updated fiel plotting script

examples/coating-flow.py

@@ -101,11 +101,14 @@ if __name__ == '__main__':
             m[0] = np.array([ms] + nvalues[:, 1].tolist(), dtype = np.complex128)
             print(x,m)
 
-            factor = 2
+            factor = 2.91*x[0][0]/x[0][-1]
+            print factor
             comment='PEC-'+basename
-            WL_units='cm'
+            WL_units=''
             #flow_total = 39
-            flow_total = 25
+            # flow_total = 23 #SV False
+            flow_total = 24
+            #flow_total = 4
             #crossplane='XZ'
             crossplane='XYZ'
             #crossplane='YZ'
@@ -115,10 +118,23 @@ if __name__ == '__main__':
             #field_to_plot='Pabs'
             #field_to_plot='Eabs'
             
-            #field_to_plot='angleEx'
-            field_to_plot='angleHy'
-            fieldplot(x[0],m[0], wl, comment, WL_units, crossplane, field_to_plot, npts,
-                      factor, flow_total, pl=0, outline_width=2.0)
+            field_to_plot='angleEx'
+            #field_to_plot='angleHy'
+
+            import matplotlib.pyplot as plt
+            fig, axs = plt.subplots(1,1)#, sharey=True, sharex=True)
+            fig.tight_layout()
+            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')
+            plt.draw()
+            plt.clf()
+            plt.close()
 
 
         print "Done!!"

examples/field-Ag-flow.py

@@ -87,5 +87,27 @@ crossplane='XZ'
 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)
 
+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(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total, is_flow_extend=False)
+
+# for ax in axs:
+#     ax.locator_params(axis='x',nbins=5)
+#     ax.locator_params(axis='y',nbins=5)
+
+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')
+
+plt.draw()
+
+#    plt.show()
+
+plt.clf()
+plt.close()

examples/field-SiAgSi-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
 #
@@ -33,43 +34,11 @@
 import scattnlay
 from scattnlay import fieldnlay
 from scattnlay import scattnlay
+from fieldplot import fieldplot
 import numpy as np
 import cmath
 
 
-def get_index(array,value):
-    idx = (np.abs(array-value)).argmin()
-    return idx
-
-#Ec = np.resize(Ec, (npts, npts)).T
-
-
-def GetFlow(scale_x, scale_z, Ec, Hc, a, b, nmax):
-    # Initial position
-    flow_x = [a]
-    flow_z = [b]
-    for n in range(0, nmax):
-        #Get the next position
-        #1. Find Poynting vector and normalize it
-        x_pos = flow_x[-1]
-        z_pos = flow_z[-1]
-        x_idx = get_index(scale_x, x_pos)
-        z_idx = get_index(scale_z, z_pos)
-        S=np.cross(Ec[npts*z_idx+x_idx], np.conjugate(Hc[npts*z_idx+x_idx]) ).real
-        Snorm=S/np.linalg.norm(S)
-        #2. Evaluate displacement = half of the discrete and new position
-        dpos = abs(scale_z[0]-scale_z[1])/4.0
-        dx = dpos*Snorm[0]
-        dz = dpos*Snorm[2]
-        x_pos = x_pos+dx
-        z_pos = z_pos+dz
-        #3. Save result
-        flow_x.append(x_pos)
-        flow_z.append(z_pos)
-    return flow_x, flow_z
-
-
-
 epsilon_Si = 13.64 + 0.047j
 epsilon_Ag = -28.05 + 1.525j
 
@@ -103,135 +72,56 @@ outer_r = inner_r+outer_width
 # n2 = 0.565838 + 7.23262j
 nm = 1.0
 
-x = np.ones((1, 3), dtype = np.float64)
-x[0, 0] = 2.0*np.pi*core_r/WL
-x[0, 1] = 2.0*np.pi*inner_r/WL
-x[0, 2] = 2.0*np.pi*outer_r/WL
+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((1, 3), dtype = np.complex128)
-m[0, 0] = index_Si/nm
-m[0, 1] = index_Ag/nm
-m[0, 2] = index_Si/nm
+m = np.ones((3), dtype = np.complex128)
+m[0] = index_Si/nm
+m[1] = index_Ag/nm
+m[2] = index_Si/nm
 
 print "x =", x
 print "m =", m
 
-npts = 241
-
+npts = 501
 factor=2.2
-scan = np.linspace(-factor*x[0, 2], factor*x[0, 2], npts)
-
-coordX, coordZ = np.meshgrid(scan, scan)
-coordX.resize(npts*npts)
-coordZ.resize(npts*npts)
-coordY = np.zeros(npts*npts, dtype = np.float64)
-
-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)
-Er = np.absolute(E)
-Hr = np.absolute(H)
-
-# |E|/|Eo|
-Eabs = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
-Ec = E[0, :, :]
-Hc = H[0, :, :]
-Eangle = np.angle(E[0, :, 0])/np.pi*180
-
-P=[]
-for n in range(0, len(E[0])):
-    P.append(np.linalg.norm( np.cross(Ec[n], np.conjugate(Hc[n]) ).real/2 ))
-
-Habs= np.sqrt(Hr[0, :, 0]**2 + Hr[0, :, 1]**2 + Hr[0, :, 2]**2)
-Hangle = np.angle(H[0, :, 1])/np.pi*180
-
-
-
-try:
-    import matplotlib.pyplot as plt
-    from matplotlib import cm
-    from matplotlib.colors import LogNorm
-
-    min_tick = 0.0
-    max_tick = 1.0
-
-    # Eabs_data = np.resize(P, (npts, npts)).T
-    Eabs_data = np.resize(Eabs, (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)#, sharey=True, sharex=True)
-    #fig.tight_layout()
-    # Rescale to better show the axes
-    scale_x = np.linspace(min(coordX)*WL/2.0/np.pi/nm, max(coordX)*WL/2.0/np.pi/nm, npts)
-    scale_z = np.linspace(min(coordZ)*WL/2.0/np.pi/nm, max(coordZ)*WL/2.0/np.pi/nm, npts)
-
-    # Define scale ticks
-    min_tick = min(min_tick, np.amin(Eabs_data))
-    max_tick = max(max_tick, np.amax(Eabs_data))
-    #max_tick = 5
-    # scale_ticks = np.power(10.0, np.linspace(np.log10(min_tick), np.log10(max_tick), 6))
-    scale_ticks = np.linspace(min_tick, max_tick, 11)
-
-    # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
-    ax.set_title('Eabs')
-    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)
-
-    plt.xlabel('Z, nm')
-    plt.ylabel('X, nm')
-
-    # This part draws the nanoshell
-    from matplotlib import patches
-    s1 = patches.Arc((0, 0), 2.0*core_r, 2.0*core_r,  angle=0.0, zorder=2,
-                     theta1=0.0, theta2=360.0, linewidth=1, color='black')
-    s2 = patches.Arc((0, 0), 2.0*inner_r, 2.0*inner_r, angle=0.0, zorder=2,
-                     theta1=0.0, theta2=360.0, linewidth=1, color='black')
-    s3 = patches.Arc((0, 0), 2.0*outer_r, 2.0*outer_r, angle=0.0, zorder=2,
-                     theta1=0.0, theta2=360.0, linewidth=1, color='black')
-    ax.add_patch(s1)
-    ax.add_patch(s2) 
-    ax.add_patch(s3) 
-
-    from matplotlib.path import Path
-    #import matplotlib.patches as patches
-
-    flow_total = 21
-    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), npts*12)
-        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='white')
-        ax.add_patch(patch)
-
- 
-    plt.savefig("SiAgSi-flow.png")
-    plt.draw()
-
-    plt.show()
-
-    plt.clf()
-    plt.close()
-finally:
-    print("Qabs = "+str(Qabs));
-#
+flow_total = 21
+
+
+crossplane='XZ'
+#crossplane='YZ'
+#crossplane='XY'
+
+# Options to plot: Eabs, Habs, Pabs, angleEx, angleHy
+field_to_plot='Eabs'
+#field_to_plot='angleEx'
+comment='SiAgSi-absorber-flow'
+WL_units='nm'
+
+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)
+
+#fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total, is_flow_extend=False)
+
+# for ax in axs:
+#     ax.locator_params(axis='x',nbins=5)
+#     ax.locator_params(axis='y',nbins=5)
+
+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')
+
+plt.draw()
+
+#    plt.show()
+
+plt.clf()
+plt.close()
 
 

examples/fieldplot.py

@@ -39,6 +39,7 @@ 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'::
 
@@ -59,66 +60,70 @@ def angle_between(v1, v2):
             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 = x[-1] / 3
+    min_step = x[0] / 2000
 #    max_step = min_step
-    step = min_step*2.0
+    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)
+    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)
+    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
+        if (count > 4000):  # Limit length of the absorbed power streamlines
             break
-        if step<max_step:
-                step = step*2.0
+        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
+            # 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)
+            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)):
+            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
+                    Ec[i] = 0 + 0j
                 if abs(Hc[i]) < Hth:
-                    Hc[i] = 0+0j
+                    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:
+            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
+            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];
+        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)
@@ -129,216 +134,201 @@ def GetFlow3D(x0, y0, z0, max_length, max_angle, x, m, pl):
 ###############################################################################
 def GetField(crossplane, npts, factor, x, m, pl):
     """
-    crossplane: XZ, YZ, XY, or XYZ (half is XZ, half is YZ)
+    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':
+    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
-    elif crossplane=='YZ':
+        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
-    elif crossplane=='XY':
+        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
-    elif crossplane=='XYZ':
-        coordX, coordZ = np.meshgrid(scan, scan)
-        coordY, coordZ = np.meshgrid(scan, scan)
-        half=npts//2
-        # coordX = np.copy(coordX)
-        # coordY = np.copy(coordY)
-        coordX[:,:half]=0
-        coordY[:,half:]=0
-        coordX.resize(npts*npts)
-        coordY.resize(npts*npts)
-        coordZ.resize(npts*npts)
-        
+        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])), range(0, len(E[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
+    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 = GetField(crossplane, npts, factor, x, m, pl)
-    scan = np.linspace(-factor*x[-1], factor*x[-1], npts)
-    coordX1, coordZ1 = np.meshgrid(scan, scan)
 
+
+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=' '):
+    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^*)$'
-            label = r'$\left |E \times H\right|$'
+            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 
+            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 
+            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 
+            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 
+            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(coordX1)*WL/2.0/np.pi, max(coordX1)*WL/2.0/np.pi, npts)
-        # scale_z = np.linspace(min(coordZ1)*WL/2.0/np.pi, max(coordZ1)*WL/2.0/np.pi, npts)
-        scale_x = np.linspace(-factor*x[-1]*WL/2.0/np.pi, factor*x[-1]*WL/2.0/np.pi, npts)
-        scale_z = np.linspace(-factor*x[-1]*WL/2.0/np.pi, factor*x[-1]*WL/2.0/np.pi, npts)
+        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)])
+        #min_tick = 0.1
         max_tick = np.amax(Eabs_data[~np.isnan(Eabs_data)])
-        scale_ticks = np.linspace(min_tick, max_tick, 6)
-
+        #max_tick = 60
+        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)
-        my_cmap = cm.jet
-        if not (field_to_plot == 'angleEx' or field_to_plot == 'angleHy'):
-            my_cmap.set_under('w')
-        cax = ax.imshow(Eabs_data, interpolation = 'nearest', cmap = my_cmap,
-                        origin = 'lower'
-                        , vmin = min_tick+max_tick*1e-15, vmax = max_tick
-                        , extent = (min(scale_x), max(scale_x), min(scale_z), max(scale_z))
-                        #,norm = LogNorm()
+        # 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))
+                        # ,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)
+        cbar = fig.colorbar(cax, ticks=[a for a in scale_ticks], ax=ax)
+        # vertically oriented colorbar
+        cbar.ax.set_yticklabels(['%3.2f' % (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)
-        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=='XYZ':
-            plt.xlabel('Z, '+WL_units)
-            plt.ylabel('Y:X, '+WL_units)
-        elif crossplane=='XY':
-            plt.xlabel('Y, '+WL_units)
-            plt.ylabel('X, '+WL_units)
-
-
+        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)
+        elif crossplane == 'YZ':
+            ax.set_xlabel('Z, ' + WL_units, labelpad=lp1)
+            ax.set_ylabel('Y, ' + WL_units, labelpad=lp2)
+        elif crossplane == 'XY':
+            ax.set_xlabel('Y, ' + WL_units, labelpad=lp1)
+            ax.set_ylabel('X, ' + WL_units, labelpad=lp2)
         # # This part draws the nanoshell
         from matplotlib import patches
         from matplotlib.path import Path
-        x_edge = (x[-1], x[0])
-        for xx in x_edge:
-            r= xx*WL/2.0/np.pi
-            s1 = patches.Arc((0, 0), 2.0*r, 2.0*r,  angle=0.0, zorder=1.8,
+        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)
-        if (not crossplane=='XY') and flow_total>0:
+        #
+        # 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, f = 0, 0, 0, 0
-            max_length=factor*x[-1]*8
-            #max_length=factor*x[-1]*4
-            max_angle = np.pi/160
+            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*2+1)
+                rg = range(0, flow_total * 5 + 1)
             else:
-                rg = range(0,flow_total)
+                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
+                if crossplane == 'XZ':
+                    if is_flow_extend:
+                        x0 = min_SP * 2 + flow * step_SP
                     else:
-                        y0 = f
-                z0 = min_SP
-
-                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
+                        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:
-                        flow_z_plot = flow_zSP*WL/2.0/np.pi
-                        flow_f_plot = flow_ySP*WL/2.0/np.pi
+                        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)
+                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)
+                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')
-        bbox_props = dict(boxstyle="round,pad=0.3", fc="w", ec="w", lw=2)
-        if crossplane=='XYZ':
-            ax.annotate('E-k', xy=(0.96, 0.96), xycoords='axes fraction', fontsize=16,
-                horizontalalignment='right', verticalalignment='top',
-                bbox=bbox_props)
-            ax.annotate('H-k', xy=(0.96, 0.04), xycoords='axes fraction', fontsize=16,
-                horizontalalignment='right', verticalalignment='bottom',
-                bbox=bbox_props)
-            ax.axhline(y=0.0, ls='--', dashes=[5,3], color='gray', lw=1.5)
+#                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+".png")
-                    +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));
+        terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(
+            np.array([x]), np.array([m]))
+        print("Qabs = " + str(Qabs))
     #
-
-