skew debug poynting

DiffractionGratings/grating3D.pro

@@ -193,6 +193,7 @@ Function{
   
   // Bloch phase shifts
   skx1[] =  k1x[];
+  // sky1[] = -k0*n1[]*Sin[theta0]*Sin[phi0+xsi];
   sky1[] = -k0*n1[]*Sin[theta0]*Sin[phi0+xsi];
 
   dephX[] = Exp[I[]*skx1[]*period_x];
@@ -201,17 +202,16 @@ Function{
 	// Fourier coefficients variables
   Nb_ordre = 2*Nmax+1;
   For i In {0:Nb_ordre-1}
-    alpha~{i}[]      = -k1x[] + 2*Pi/period_x * (i-Nmax);
-    expialphax~{i}[] = Exp[I[]*alpha~{i}[]*X[]];
-  EndFor
     For j In {0:Nb_ordre-1}
-    beta~{j}[]      = -k1y[]  + 2*Pi/period_y * (j-Nmax);
-    expibetay~{j}[] = Exp[I[]*beta~{j}[]*Y[]];
+      alpha~{i}~{j}[] = -k1x[] + 2*Pi/period_x*(i-Nmax);
+      beta~{i}~{j}[]  = -k1y[] + 2*Pi/period_y*(j-Nmax)/Cos[xsi] - 2*Pi/period_x*(i-Nmax)*Tan[xsi];
+      expialphaxy~{i}~{j}[] = Exp[I[]*(alpha~{i}~{j}[]*X[]+beta~{i}~{j}[]*Y[])];
+    EndFor
   EndFor
   For i In {0:Nb_ordre-1}
     For j In {0:Nb_ordre-1}
-      gammar~{i}~{j}[] = Sqrt[k0^2*epsr1[] - alpha~{i}[]^2 - beta~{j}[]^2];
-      gammat~{i}~{j}[] = Sqrt[k0^2*epsr2[] - alpha~{i}[]^2 - beta~{j}[]^2];
+      gammar~{i}~{j}[] = Sqrt[k0^2*epsr1[] - alpha~{i}~{j}[]^2 - beta~{i}~{j}[]^2];
+      gammat~{i}~{j}[] = Sqrt[k0^2*epsr2[] - alpha~{i}~{j}[]^2 - beta~{i}~{j}[]^2];
     EndFor
   EndFor
 
@@ -342,6 +342,8 @@ PostProcessing {
       { Name E1perx ; Value { Local { [ CompX[E1[]*Exp[-I[]*(skx1[]*X[]+sky1[]*Y[])]] ]; In Omega; Jacobian JVol; } } }
       { Name E1pery ; Value { Local { [ CompY[E1[]*Exp[-I[]*(skx1[]*X[]+sky1[]*Y[])]] ]; In Omega; Jacobian JVol; } } }
       
+      { Name Etotpery ; Value { Local { [ CompY[({u}+E1[])*Exp[-I[]*(k1x[]*X[]+k1y[]*Y[])]] ]; In Omega; Jacobian JVol; } } }
+      
       { Name Etot   ; Value { Local { [ {u}+E1[]  ]; In Omega; Jacobian JVol; } } }
       { Name Htot   ; Value { Local { [ H1[]-I[]/(mur[]*mu0*om0)*{Curl u}]; In Omega; Jacobian JVol; } } }
       { Name Htotx  ; Value { Local { [CompX[ H1[]-I[]/(mur[]*mu0*om0)*{Curl u}]]; In Omega; Jacobian JVol; } } }
@@ -365,28 +367,28 @@ PostProcessing {
 
       For i In {0:Nb_ordre-1}
         For j In {0:Nb_ordre-1}
-          { Name int_x_t~{i}~{j} ; Value { Integral { [ CompX[{u}+E1[] ]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
-          { Name int_y_t~{i}~{j} ; Value { Integral { [ CompY[{u}+E1[] ]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
-          { Name int_x_r~{i}~{j} ; Value { Integral { [ CompX[{u}+E1d[]]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
-          { Name int_y_r~{i}~{j} ; Value { Integral { [ CompY[{u}+E1d[]]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
+          { Name int_x_t~{i}~{j} ; Value { Integral { [ CompX[{u}+E1[] ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
+          { Name int_y_t~{i}~{j} ; Value { Integral { [ CompY[{u}+E1[] ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
+          { Name int_x_r~{i}~{j} ; Value { Integral { [ CompX[{u}+E1d[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
+          { Name int_y_r~{i}~{j} ; Value { Integral { [ CompY[{u}+E1d[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
           { Name eff_t~{i}~{j}   ; Value { Term{Type Global; [
-                 1/(gammat~{i}~{j}[]*-k1z[]) *((gammat~{i}~{j}[]^2+alpha~{i}[]^2)*SquNorm[$int_x_t~{i}~{j}]+
-                                                       (gammat~{i}~{j}[]^2+ beta~{j}[]^2)*SquNorm[$int_y_t~{i}~{j}]+
-                                                       2*alpha~{i}[]*beta~{j}[]*Re[$int_x_t~{i}~{j}*Conj[$int_y_t~{i}~{j}]]) ] ; In SurfIntBot ; } } }
+                 1/(gammat~{i}~{j}[]*-k1z[]) * ( (gammat~{i}~{j}[]^2+alpha~{i}~{j}[]^2)*SquNorm[$int_x_t~{i}~{j}]+
+                                              (gammat~{i}~{j}[]^2+  beta~{j}~{j}[]^2)*SquNorm[$int_y_t~{i}~{j}]+
+                                              2*alpha~{i}~{j}[]*beta~{i}~{j}[]*Re[$int_x_t~{i}~{j}*Conj[$int_y_t~{i}~{j}]] ) ] ; In SurfIntBot ; } } }
           { Name eff_r~{i}~{j}   ; Value { Term{Type Global; [
-                 1/(gammar~{i}~{j}[]*-k1z[])*((gammar~{i}~{j}[]^2+alpha~{i}[]^2)*SquNorm[$int_x_r~{i}~{j}]+
-                                                      (gammar~{i}~{j}[]^2+ beta~{j}[]^2)*SquNorm[$int_y_r~{i}~{j}]+
-                                                      2*alpha~{i}[]*beta~{j}[]*Re[$int_x_r~{i}~{j}*Conj[$int_y_r~{i}~{j}]]) ] ; In SurfIntTop ; } } }
+                 1/(gammar~{i}~{j}[]*-k1z[])* ( (gammat~{i}~{j}[]^2+alpha~{i}~{j}[]^2)*SquNorm[$int_x_r~{i}~{j}]+
+                                                      (gammar~{i}~{j}[]^2+ beta~{i}~{j}[]^2)*SquNorm[$int_y_r~{i}~{j}]+
+                                                      2*alpha~{i}~{j}[]*beta~{i}~{j}[]*Re[$int_x_r~{i}~{j}*Conj[$int_y_r~{i}~{j}]]) ] ; In SurfIntTop ; } } }
           { Name numbering_ij~{i}~{j}   ; Value { Term{Type Global; [Vector[i-Nmax,j-Nmax,0]] ; In SurfIntBot ; } } }
         EndFor
       EndFor
 
-      For i In {0:Nb_ordre-1}
-        { Name alpha~{i}       ; Value { Term{Type Global; [alpha~{i}[]] ; In SurfIntBot ; } } }
-      EndFor
-      For j In {0:Nb_ordre-1}
-        { Name beta~{j}   ; Value { Term{Type Global; [beta~{j}[]] ; In SurfIntBot ; } } }
-      EndFor
+      // For i In {0:Nb_ordre-1}
+      //   For j In {0:Nb_ordre-1}
+      //     { Name alpha~{i}~{j} ; Value { Term{Type Global; [alpha~{i}[]-alpha~{i}[]] ; In SurfIntBot ; } } }
+      //     { Name beta~{j}~{j}  ; Value { Term{Type Global; [beta~{j}[]] ; In SurfIntBot ; } } }
+      //   EndFor
+      // EndFor
     }
   }
 }
@@ -443,8 +445,9 @@ PostOperation {
       // Print [ upery  , OnElementsOf Omega, File StrCat[myDir,"upery.pos"]];
       // Print [ uperz  , OnElementsOf Omega, File StrCat[myDir,"uperz.pos"]];
 
-      Print [ E1perx , OnElementsOf Omega, File StrCat[myDir,"E1perx.pos"]];
-      Print [ E1pery , OnElementsOf Omega, File StrCat[myDir,"E1pery.pos"]];
+      // Print [ E1perx , OnElementsOf Omega, File StrCat[myDir,"E1perx.pos"]];
+      // Print [ E1pery , OnElementsOf Omega, File StrCat[myDir,"E1pery.pos"]];
+      Print [ Etotpery , OnElementsOf Omega, File StrCat[myDir,"Etotpery.pos"]];
       Print [ Poy_tot , OnPlane { {0.5*(-period_x-dys), -dyc/2,(hh_L_6+hh_L_5)/2}
                                   {0.5*( period_x-dys), -dyc/2,(hh_L_6+hh_L_5)/2}
                                   {0.5*(-period_x+dys),  dyc/2,(hh_L_6+hh_L_5)/2} }

DiffractionGratings/grating3D_postplot.py

@@ -10,7 +10,7 @@ def dtrap_poy(fname_in,nx,ny):
     return np.trapz(temp,y_export2D[0,:]) #[x_export2D,y_export2D,poy_y_grid_re] #
 
 myDir = sys.argv[1]
-
+fact = np.cos(float(sys.argv[2])*np.pi/180)**2
 intpoyz_tot =  dtrap_poy(myDir+'/Poy_tot_gd.pos',50,50) #*np.cos(float(sys.argv[2])*np.pi/180)
 intpoyz_ref =  dtrap_poy(myDir+'/Poy_ref_gd.pos',50,50) #*np.cos(float(sys.argv[2])*np.pi/180)
 intpoyz_inc =  dtrap_poy(myDir+'/Poy_inc_gd.pos',50,50) #*np.cos(float(sys.argv[2])*np.pi/180)
@@ -24,6 +24,8 @@ Ascat2      = np.loadtxt(myDir+'/temp-Q_scat2.txt')[1]
 
 Rnm = np.loadtxt(myDir+'/eff_r.txt',ndmin=2)[:,1] + 1j*np.loadtxt(myDir+'/eff_r.txt',ndmin=2)[:,2]
 Tnm = np.loadtxt(myDir+'/eff_t.txt',ndmin=2)[:,1] + 1j*np.loadtxt(myDir+'/eff_t.txt',ndmin=2)[:,2]
+Rnm/=fact
+Tnm/=fact
 Q = [np.loadtxt(myDir+'/temp-Q_L_%g.txt'%k,ndmin=2)[:,1] for k in range(2,7)]
 Q.append(np.loadtxt(myDir+'/temp-Q_scat.txt',ndmin=2)[:,1])
 Q=np.array(Q)