Add cloak in glass example, coating flow script update

examples/coating-flow.py

@@ -99,14 +99,15 @@ if __name__ == '__main__':
 
             x[0] = 2.0*np.pi*np.array(r, dtype = np.float64)/wl
             m[0] = np.array([ms] + nvalues[:, 1].tolist(), dtype = np.complex128)
+            print(x,m)
 
             factor = 2
             comment='PEC-'+basename
             WL_units='cm'
-            flow_total = 39
-            #flow_total = 5
-            #crossplane='XZ'
-            crossplane='YZ'
+            #flow_total = 39
+            flow_total = 25
+            crossplane='XZ'
+            #crossplane='YZ'
             #crossplane='XY'
 
             # Options to plot: Eabs, Habs, Pabs, angleEx, angleHy

examples/coating-flow.sh

@@ -1,6 +1,8 @@
 #!/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.80 -f index-dv.dat -n 1501
 #python coating-flow.py -w 3.75 -t 0.80 -f index-dv.dat -n 51
 #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 1501

examples/fieldplot.py

@@ -165,7 +165,7 @@ def GetField(crossplane, npts, factor, x, m, 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]))))
+    P = np.array(map(lambda n: np.linalg.norm(np.cross(Ec[n], np.conjugate(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 ))
@@ -181,7 +181,7 @@ def fieldplot(x,m, WL, comment='', WL_units=' ', crossplane='XZ', field_to_plot=
         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'$\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 
@@ -212,7 +212,8 @@ def fieldplot(x,m, WL, comment='', WL_units=' ', crossplane='XZ', field_to_plot=
         # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
         ax.set_title(label)
         my_cmap = cm.jet
-        my_cmap.set_under('w')
+        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
@@ -293,7 +294,8 @@ def fieldplot(x,m, WL, comment='', WL_units=' ', crossplane='XZ', field_to_plot=
                 #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+".png")
+                    +field_to_plot+".pdf")
         plt.draw()
 
     #    plt.show()

examples/index-in-glass.dat

@@ -0,0 +1,6 @@
+#Layer N	Index	
+# TargetR0.3990-CoatingW00.128-FinalRCS00.4505diff-57.5%-n4-s04.817268921311-index
+1	0.7662439167558822151	
+2	2.588110333651454376	
+3	0.7397603938619315977	
+4	0.7231542770413039722