SiAgSi with powerflow lines

README.md

@@ -53,12 +53,14 @@ terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
   ```bash
 scattnlay -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] [-c comment]
   ```
+
   * Execute some of the test scripts (located in the folder 'tests/shell')
       Example:
 
   ```bash
 ./test01.sh > test01.csv
   ```
+  
 3. C++ library
 
 ```C++

tests/python/field-SiAgSi.py

@@ -34,12 +34,19 @@ import scattnlay
 from scattnlay import fieldnlay
 from scattnlay import scattnlay
 import numpy as np
+import cmath
 
 epsilon_Si = 13.64 + 0.047j
 epsilon_Ag = -28.05 + 1.525j
+
 # epsilon_Si = 2.0 + 0.047j
 # epsilon_Ag = -2.0 + 1.525j
 
+# air = 1
+# epsilon_Si = air*2
+# epsilon_Ag = air*2
+
+
 index_Si = np.sqrt(epsilon_Si)
 index_Ag = np.sqrt(epsilon_Ag)
 
@@ -75,7 +82,8 @@ print "m =", m
 
 npts = 281
 
-scan = np.linspace(-2.0*x[0, 2], 2.0*x[0, 2], npts)
+factor=4.5
+scan = np.linspace(-factor*x[0, 2], factor*x[0, 2], npts)
 
 coordX, coordZ = np.meshgrid(scan, scan)
 coordX.resize(npts*npts)
@@ -88,63 +96,116 @@ terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
 terms, E, H = fieldnlay(x, m, coord)
 print("Qabs = "+str(Qabs));
 Er = np.absolute(E)
+Hr = np.absolute(H)
 
 # |E|/|Eo|
-Eh = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
+Eabs = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
+Eangle = np.angle(E[0, :, 0])/np.pi*180
 
-result = np.vstack((coordX, coordY, coordZ, Eh)).transpose()
+Habs= np.sqrt(Hr[0, :, 0]**2 + Hr[0, :, 1]**2 + Hr[0, :, 2]**2)
+Hangle = np.angle(H[0, :, 1])/np.pi*180
+
+result = np.vstack((coordX, coordY, coordZ, Eabs)).transpose()
+result2 = np.vstack((coordX, coordY, coordZ, Eangle)).transpose()
 
 try:
     import matplotlib.pyplot as plt
     from matplotlib import cm
     from matplotlib.colors import LogNorm
 
-    min_tick = 0.0
-    max_tick = 1.0
+    # min_tick = 0.0
+    # max_tick = 1.0
 
-    edata = np.resize(Eh, (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 = plt.figure()
-    ax = fig.add_subplot(111)
+    fig, axs = plt.subplots(2,2)#, 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(edata))
-    max_tick = max(max_tick, np.amax(edata))
+    # min_tick = min(min_tick, np.amin(Eabs_data))
+    # max_tick = max(max_tick, np.amax(Eabs_data))
     # 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, 12)
+    # scale_ticks = np.linspace(min_tick, max_tick, 10)
 
     # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
-    cax = ax.imshow(edata, 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))
+    axs[0,0].set_title('Eabs')
+    cax = axs[0,0].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()
+                    )
+    axs[0,0].axis("image")
+    axs[0,1].set_title('Eangle')
+    cax = axs[0,1].imshow(Eangle_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()
+                    )
+    axs[1,0].set_title('Habs')
+    cax = axs[1,0].imshow(Habs_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()
+                    )
+    axs[1,1].set_title('Hangle')
+    cax = axs[1,1].imshow(Hangle_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()
                     )
 
     # 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)
+    # 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')
+    # 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')
-    ax.add_patch(s1)
-
-    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')
-    ax.add_patch(s2) 
-    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(s3) 
+    for m in (0,1):
+        for n in (0,1):
+            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')
+            axs[m,n].add_patch(s1)
+            axs[m,n].add_patch(s2) 
+            axs[m,n].add_patch(s3) 
+    # axs[0,0].add_patch(s1)
+    # axs[0,0].add_patch(s2)
+    # axs[0,0].add_patch(s3)
+    # axs[1,0].add_patch(s1)
+    # axs[1,0].add_patch(s2)
+    # axs[1,0].add_patch(s3)
+    # axs[0,1].add_patch(s1)
+    # axs[0,1].add_patch(s2)
+    # axs[0,1].add_patch(s3)
+    # axs[1,1].add_patch(s1)
+    # axs[1,1].add_patch(s2)
+    # axs[1,1].add_patch(s3)
+            
+    # for m in (0,1):
+    #     for n in (0,1):
+    #         print(m)
+    #         print(n)
+    #         axs[m,n].add_patch(s1)
+    #         axs[m,n].add_patch(s2) 
+    #         axs[m,n].add_patch(s3) 
     # End of drawing
     plt.savefig("SiAgSi.png")
     plt.draw()