updated HPC scripts

add_plasmon_monitors.lsf

@@ -1,6 +1,8 @@
 groupscope("::model");
 select("FDTD");
 min_mesh_step=get("min mesh step");
+z_mon = get("z");
+z_span = get("z span");
 select("sub-Au");
 sub_Au_z_max = get("z max");
 #remove old monitors
@@ -17,11 +19,24 @@ for(x=1:10) {
             addpower;
             zshift_all = z_arr(i)*1e-9;
             set("name","mon_x"+num2str(x)+"mkm_media"+num2str(media)+"_zshift"
-                +num2str(z_arr(i))+"nm");
+            +num2str(z_arr(i))+"nm");
             set("monitor type",1);  # 1 = point, 2 = linear x, 3 = linear y, 4 = linear z, 5 = 2D x-normal, 6 = 2D y-normal, 7 = 2D z-normal, 8 = 3D
             set("x",x*mkm);
             set("y",0);
             set("z",sub_Au_z_max + shift_sign*zshift_all);
         }
     }
-}
+}
+groupscope("::model");
+for(i=1:2){
+    select("monitor_"+num2str(i));
+    delete;
+    addpower;
+    set("name","monitor_"+num2str(i));
+    set("monitor type",4);  # 1 = point, 2 = linear x, 3 = linear y, 4 = linear z, 5 = 2D x-normal, 6 = 2D y-normal, 7 = 2D z-normal, 8 = 3D
+    set("x",0.15*mkm+0.650*(i-1)*mkm);
+    set("y",0);
+    set("z", z_mon);
+    set("z span", z_span);
+}
+

create_hpc_bat.lsf

@@ -5,18 +5,41 @@ slash = "/";
 if (operatingsystem == "windows") {slash = "\\";}
 
 curdir = pwd;
+if (fileexists("temppath.txt")) {
+    temppath = read("temppath.txt");
+    len = length(temppath);
+    temppath = substring(temppath,1,len-2);
+    cd(temppath);
+    if (fileexists("hostname.txt")) {
+        hostname = read("hostname.txt");
+        len = length(hostname);
+        hostname = substring(hostname,1,len-2);
+        cd(curdir);
+    } else {
+        cd(curdir);
+        break;
+    }
+} else {
+    break;
+}
 files = splitstring(dir,endl);
 isCreated = false;
-fname = curdir+slash+"exec.bat";
+fname = curdir+slash+hostname+"exec.bat";
 if (fileexists(fname)){rm(fname);}
 for(i=1:length(files)) {
     if(fileextension(files{i}) == "fsp") {
-        load(files{i});
-        status = layoutmode;
-        if (status == 1) {
-            write(fname,"\"C:\Program Files\Lumerical\FDTD\bin\fdtd-engine\" -t 32 "+files{i});
-            write(fname,"exit /b 1");
-            break;
+        lockfile =  files{i}+".lock";
+        if (!fileexists(lockfile)){
+            load(files{i});
+            status = layoutmode;
+            if (status == 1) {
+                write(fname,"IF EXIST \""+lockfile+"\" exit /b 1");
+                write(fname,"copy /y NUL "+lockfile+" >NUL");
+                write(fname,"\"C:\Program Files\Lumerical\FDTD\bin\fdtd-engine\" -t 32 "+files{i});
+                write(fname,"del "+lockfile);
+                write(fname,"exit /b 1");
+                break;
+            }
         }
     }
 }

efficiency-plasmon-plot.py

@@ -108,6 +108,9 @@ def analyze(data, dist, z_vec, wl_idx):
     kappa1= np.sqrt(k_spp**2 - eps1*k_0**2)
     kappa2= np.sqrt(k_spp**2 - eps2*k_0**2)
 
+    print(1e9*lambd[9])
+    print(1e9/kappa1[9])
+    print(1e9/kappa2[9])
     check_field_match(data_in_air, data_in_gold,wl_idx,z_vec,kappa1,kappa2,eps2)
 
     H1 = data_in_air[:,6]
@@ -166,8 +169,10 @@ file_ext="pdf"
 #dirname="Au-JC-R100-Au-JC.fsp.results"
 #dirname="Au-McPeak-R100-Si-Green.fsp.results"
 #dirname="Au-McPeak-R100-Au-McPeak.fsp.results"
-dirname="sub-Au-R100-Si-wl450-800-sep10nm.fsp.results"
-
+#dirname="sub-Au-R100-Si-wl450-800-sep10nm.fsp.results"
+#dirname="bg-Au-sub-Au-dipole-W.fsp.results"
+#dirname="bg-Au-sub-Si-dipole-W.fsp.results"
+dirname="bg-Au-sub-dipole-W.fsp.results"
 #dirname="Au-McPeak-R0.fsp.results"
 #dirname="Au-McPeak-R100-Si-Green-1500.fsp.results"
 #dirname="Au-McPeak-R100-Si-Green-1500-l.fsp.results"

hpc.bat

@@ -6,11 +6,19 @@ set "startTime=%time: =0%"
 REM path %path%;
 REM"C:\Program Files\Lumerical\FDTD\bin\fdtd-engine" -t 32 template-dipole-on-sphere-on-surf-z.fsp
 
+set host=%COMPUTERNAME%
+echo %host% > %TEMP%\hostname.txt
+
+echo exit /b 1 > %host%exec.bat
+
 echo Before loop
 :loop
+REM random timeout 10 seconds
+    set /a timeout=%RANDOM% * 11 / 32768
+    timeout %timeout%
+    echo %TEMP% > temppath.txt
     "C:\Program Files\Lumerical\FDTD\bin\fdtd-solutions" -nw -run create_hpc_bat.lsf
-    echo before run
-    call exec.bat
+    call %host%exec.bat
 if !errorlevel! gtr 0 goto loop
 
 echo after loop

main.lsf

@@ -6,6 +6,7 @@ if (status == 1) {
 } else {
 #    switchtolayout;
     save_plasmon_data;
+    mat_save;
 }
 #newproject("current");
 redrawon;

plasmon-modal-efficiency.py

@@ -9,12 +9,15 @@ import scipy.io
 c = 299792458.0
 eps_0 = 8.854187817e-12 # F/m
 pi = np.pi
-verbose = 6
 # r of monitor
 #r = 146.513e-9
 r = 800e-9
-#debug = True
+# debug = True
+# verbose = 7
+
 debug = False
+verbose = 6
+
 
 def read_data_mat(fname):
     #data = "z, dip.power, Ex, Ey, Ez, Hx, Hy, Hz, n_Au"
@@ -23,6 +26,7 @@ def read_data_mat(fname):
     
     lambd = np.reshape(mat["lambda"],(-1))
     dippower = np.reshape(mat["dippower"],(-1))
+    vacpower = np.reshape(mat["vacpower"],(-1))
     z = np.reshape(mat["z"],(-1))
     n_Au = np.reshape(mat["n_fdtd"],(-1))
     onez = np.ones((len(z)))
@@ -38,13 +42,13 @@ def read_data_mat(fname):
         fdata = np.vstack(( z.astype(np.complex128), dippower[i]*onez.astype(np.complex128)
                             ,E1[:,0,i], E1[:,1,i], E1[:,2,i]
                             ,H1[:,0,i], H1[:,1,i], H1[:,2,i]
-                            ,n_Au[i]*onez.astype(np.complex128)
+                            ,n_Au[i]*onez.astype(np.complex128), dippower[i]*onez.astype(np.complex128)
         ))
         data1[lambd[i]]=fdata
         fdata = np.vstack(( z.astype(np.complex128), dippower[i]*onez.astype(np.complex128)
                             ,E2[:,0,i], E2[:,1,i], E2[:,2,i]
                             ,H2[:,0,i], H2[:,1,i], H2[:,2,i]
-                            ,n_Au[i]*onez.astype(np.complex128)
+                            ,n_Au[i]*onez.astype(np.complex128), dippower[i]*onez.astype(np.complex128)
         ))
         data2[lambd[i]]=fdata
         if debug: break
@@ -107,18 +111,24 @@ def get_WLs_idx(WLs, data):
 #             print("H0 air  (%5.4g %+5.4gj)"%(np.real(H1_0), np.imag(H1_0)),
 #                   " from H1 (%5.4g %+5.4gj)"%(np.real(H1[i]), np.imag(H1[i])))
 
+def cross(a, b):
+    c = [a[1]*b[2] - a[2]*b[1],
+         a[2]*b[0] - a[0]*b[2],
+         a[0]*b[1] - a[1]*b[0]]
+
+    return c
 
 def analyze(data,wl):
     # print(data[0,:]) # all z values
-    #data = "z, dip.power, Ex, Ey, Ez, Hx, Hy, Hz, n_Au"
-    #        0, 1        , 2 , 3 , 4 , 5 , 6 , 7 , 8   "
+    #data = "z, dip.power, Ex, Ey, Ez, Hx, Hy, Hz, n_Au, vac. power"
+    #        0, 1        , 2 , 3 , 4 , 5 , 6 , 7 , 8   , 9"
     lambd = wl
     omega = 2*pi*c/lambd
-
     eps_d = complex(1) # air, z>0
     eps_m = data[8,0]**2 # metal, z<0
 
     dip_power = data[1,0]
+    vac_power = data[9,0]
 
     z = data[0,:]
 
@@ -135,8 +145,6 @@ def analyze(data,wl):
         raise
     
 
-    Ex = data[2,:]
-
     Ex_m = data[2,idx_m][0]
     Ey_m = data[3,idx_m][0]
     Ez_m = data[4,idx_m][0]
@@ -145,7 +153,7 @@ def analyze(data,wl):
     Hz_m = data[7,idx_m][0]
     E_m = np.transpose(np.array([Ex_m,Ey_m,Ez_m]))
     H_m = np.transpose(np.array([Hx_m,Hy_m,Hz_m]))
-
+    
     Ex_d = data[2,idx_d][0]
     Ey_d = data[3,idx_d][0]
     Ez_d = data[4,idx_d][0]
@@ -161,13 +169,31 @@ def analyze(data,wl):
     chi_d = np.sqrt( eps_d*k_0**2 - k_sp**2 ) # desc. after eq6c, supmat
     chi_m = np.sqrt( eps_m*k_0**2 - k_sp**2 ) # desc. after eq6c, supmat
 
-    h_sp_d = np.exp(1j*chi_d*z_d) # eq6a, supmat
-    e_sp_x_d = chi_d/(omega*eps_0*eps_d)*np.exp(1j*chi_d*z_d) # eq6b, supmat
-    e_sp_z_d = k_sp/(omega*eps_0*eps_d)*np.exp(1j*chi_d*z_d) # eq6c, supmat
+    # # TODO!!! alt
+    # chi_d = np.sqrt( k_sp**2 - eps_d*k_0**2 ) # desc. after eq6c, supmat
+    # chi_m = np.sqrt( k_sp**2 - eps_m*k_0**2 ) # desc. after eq6c, supmat
+
+
+    mul = -1 #TODO!!! minus???
+    h_sp_d = np.exp(mul*1j*chi_d*z_d) # eq6a, supmat 
+    e_sp_x_d = chi_d/(omega*eps_0*eps_d)*np.exp(mul*1j*chi_d*z_d) # eq6b, supmat
+    e_sp_z_d = -k_sp/(omega*eps_0*eps_d)*np.exp(mul*1j*chi_d*z_d) # eq6c, supmat
 
     h_sp_m = np.exp(1j*-chi_m*z_m) # eq6a, supmat
     e_sp_x_m = -chi_m/(omega*eps_0*eps_m)*np.exp(1j*-chi_m*z_m) # eq6b, supmat
-    e_sp_z_m = k_sp/(omega*eps_0*eps_m)*np.exp(1j*-chi_m*z_m) # eq6c, supmat
+    e_sp_z_m = -k_sp/(omega*eps_0*eps_m)*np.exp(1j*-chi_m*z_m) # eq6c, supmat
+
+    if verbose > 8:
+        print("WL =",1e9*wl)
+        # if 1e9*wl > 763 and 1e9*wl < 766:
+        print(1j*1e9/chi_d)
+        print(1j*1e9/chi_m)
+        print("_d")
+        print(1e9*z_d)
+        print(h_sp_d)
+        print("_m")
+        print(1e9*z_m)
+        print(h_sp_m)
 
     zero_m = np.zeros(len(h_sp_m))
     zero_d = np.zeros(len(h_sp_d))
@@ -188,14 +214,22 @@ def analyze(data,wl):
 
     # E_m H_m E_d H_d
 
-    N_sp_0 = (((-1)**0) * (4.0j/(omega*eps_0*k_sp))
+    N_sp_0 = ( ((-1)**0) * (4.0j / (omega*eps_0*k_sp))
                 * (eps_d**2 - eps_m**2) / ((eps_m*eps_d)**(3/2)) )
     
     tmp_m = np.cross(E_minus_sp_0_m,H_m) - np.cross(E_m, H_minus_sp_0_m)
+    #tmp_m = np.cross(E_minus_sp_0_m,np.conj(H_m)) + np.cross(np.conj(E_m), H_minus_sp_0_m)
     radail_pojeciton_m = np.transpose(tmp_m)[0]
     integrand_m = (2*pi/N_sp_0)*radail_pojeciton_m*r 
+    # a = E_minus_sp_0_m[0]*1e9
+    # print("\n\nE_minus:", a )
+    # b = H_m[0]*1e9
+    # print("H_m:",b)
+    # print("\ncross:",np.cross(E_minus_sp_0_m*1e9,H_m*1e9)[0])
+    # print("cross:",np.array(cross(a,b)))
 
     tmp_d = np.cross(E_minus_sp_0_d,H_d) - np.cross(E_d, H_minus_sp_0_d)
+    #tmp_d = np.cross(E_minus_sp_0_d,np.conj(H_d)) + np.cross(np.conj(E_d), H_minus_sp_0_d)
     radail_pojeciton_d = np.transpose(tmp_d)[0]
     integrand_d = (2*pi/N_sp_0)*radail_pojeciton_d*r 
     
@@ -243,13 +277,15 @@ def analyze(data,wl):
 file_ext="pdf"
 
 
-def main ():
+def main (monitor_index):
     #dirname="bigourdan-Au-sub-dipole-W.fsp.1D.monitor_1.results"
-    dirname="bigourdan-Au-sub-dipole-W-2mon.fsp.1D.monitor_2.results"
-    data = read_data(dirname)
-    filename = 'bigourdan-Au-sub-dipole-W-2mon.fsp.1D.mat'
+    #dirname="bigourdan-Au-sub-dipole-W-2mon.fsp.1D.monitor_2.results"
+    #data = read_data(dirname)
+    #filename = 'bigourdan-Au-sub-dipole-W-2mon.fsp.1D.mat'
+    filename="bg-Au-sub-dipole-W.fsp.1D.mat"
+
     data2 = read_data_mat(filename)
-    r,data = data2[0]
+    r,data = data2[monitor_index]
     if verbose > 5:
         print("r =",r)
     WLs = []
@@ -257,16 +293,17 @@ def main ():
     for wl in data:
         WLs.append(wl)
         A2.append(analyze(data[wl],wl))
-
-        #print(WLs)
     WLs1 = np.array(WLs)
     A21 = np.array(A2)
+    #return
+
     # #dirname="bigourdan-Au-sub-Cyl-dipole-W.fsp.1D.monitor_1.results"
     # dirname="bigourdan-Au-sub-Cyl-dipole-W-2mon.fsp.1D.monitor_2.results"
     # data = read_data(dirname)
-    filename = 'bigourdan-Au-sub-Cyl-dipole-W-2mon.fsp.1D.mat'
+    filename = 'bg-Au-sub-Au-dipole-W.fsp.1D.mat'
+    #filename = 'bg-Au-sub-Si-dipole-W.fsp.1D.mat'
     data2 = read_data_mat(filename)
-    r,data = data2[0]
+    r,data = data2[monitor_index]
     WLs = []
     A2 = []
     for wl in data:
@@ -280,15 +317,21 @@ def main ():
     # data = np.vstack((WLs,A2))
     # print(np.sort(data))
     
-    plt.plot(WLs1*1e9, A21*275, color="black",label="x 275, no ant.")
-    plt.plot(WLs2*1e9, A22, color="red", label="with antena")
+    plt.plot(WLs1*1e9, A22/A21, color="black",label="eff.")
+    # plt.plot(WLs1*1e9, A21, color="black",label="x 275, no ant.")
+    # plt.plot(WLs2*1e9, A22, color="red", label="with antena")
     plt.legend()
     plt.xlabel(r'$\lambda$, nm')
-    plt.ylim(0,0.5)
-    plt.ylabel(r'$|A_{sp}|^2$',labelpad=-1)
-    #plt.title(dirname)
-    plt.savefig(dirname+"_A2."+file_ext)
+    # if monitor_index == 0:
+    #     plt.ylim(0,2)
+    # else:
+    #     plt.ylim(0,0.32)
+    #plt.ylim(0,2)
+    plt.ylabel(r'$^{|A^{ant}_{sp}|^2}/_{|A^0_{sp}|^2}$',labelpad=1,size=14)
+    plt.title("r = %4.1f nm"%(r*1e9))
+    plt.savefig(filename+str(monitor_index)+"_A2."+file_ext)
     plt.clf()
     plt.close()
 
-main()
+main(0)
+main(1)