skew debug poynting

DiffractionGratings/grating3D.pro

@@ -139,19 +139,21 @@ Function{
 	epsr_annex[L_1]          = epsr[];
 	epsr_annex[L_6]          = epsr[];
 
-  //// Reference Field solution of annex problem (simple diopter)  
-  k1[]    = k0*n1[]*Vector[-Sin[theta0]*Cos[phi0], -Sin[theta0]*Sin[phi0], -Cos[theta0]];
-  k2[] = Vector[
-            CompX[k1[]],
-            CompY[k1[]],
-           -Sqrt[k0^2*epsr2[]-CompX[k1[]]^2-CompY[k1[]]^2]
-        ];
-  k1r[]  = Vector[CompX[k1[]], CompY[k1[]], -CompZ[k1[]]];
-
-  rs[] = (CompZ[k1[]]-CompZ[k2[]])/(CompZ[k1[]]+CompZ[k2[]]);
-  ts[] =            2.*CompZ[k1[]] /(CompZ[k1[]]+CompZ[k2[]]);
-  rp[] = (CompZ[k1[]]*epsr2[]-CompZ[k2[]]*epsr1[])/(CompZ[k1[]]*epsr2[]+CompZ[k2[]]*epsr1[]);
-  tp[] =                   (2.*CompZ[k1[]]*epsr2[])/(CompZ[k1[]]*epsr2[]+CompZ[k2[]]*epsr1[]);
+  //// Reference Field solution of annex problem (simple diopter)
+  k1x[] = -k0*n1[]*Sin[theta0]*Cos[phi0];
+  k1y[] = -k0*n1[]*Sin[theta0]*Sin[phi0];
+  k1z[] = -k0*n1[]*Cos[theta0];
+  k2x[] =  k1x[];
+  k2y[] =  k1y[];
+  k2z[] = -Sqrt[k0^2*epsr2[]-k1x[]^2-k1y[]^2];
+  k1[]  = Vector[k1x[],k1y[], k1z[]];
+  k2[]  = Vector[k2x[],k2y[], k2z[]];
+  k1r[] = Vector[k1x[],k1y[],-k1z[]];
+
+  rs[] = (k1z[]-k2z[])/(k1z[]+k2z[]);
+  ts[] =    2.*k1z[] /(k1z[]+k2z[]);
+  rp[] = (k1z[]*epsr2[]-k2z[]*epsr1[])/(k1z[]*epsr2[]+k2z[]*epsr1[]);
+  tp[] =            (2.*k1z[]*epsr2[])/(k1z[]*epsr2[]+k2z[]*epsr1[]);
 
   spol[]    = Vector[Sin[phi0],-Cos[phi0],0];
   AmplEis[] =      spol[];
@@ -190,17 +192,20 @@ Function{
   source[] = (om0/cel)^2*(epsr[]-epsr_annex[])*E1[];
   
   // Bloch phase shifts
-  dephX[] = Exp[I[]*CompX[k1[]]*period_x];
-	dephY[] = Exp[I[]*CompY[k1[]]*period_y];
+  skx1[] =  k1x[];
+  sky1[] = -k0*n1[]*Sin[theta0]*Sin[phi0+xsi];
+
+  dephX[] = Exp[I[]*skx1[]*period_x];
+	dephY[] = Exp[I[]*sky1[]*period_y];
 
 	// Fourier coefficients variables
   Nb_ordre = 2*Nmax+1;
   For i In {0:Nb_ordre-1}
-    alpha~{i}[]      = -CompX[k1[]] + 2*Pi/period_x * (i-Nmax);
+    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}[]      = -CompY[k1[]]  + 2*Pi/period_y * (j-Nmax);
+    beta~{j}[]      = -k1y[]  + 2*Pi/period_y * (j-Nmax);
     expibetay~{j}[] = Exp[I[]*beta~{j}[]*Y[]];
   EndFor
   For i In {0:Nb_ordre-1}
@@ -227,7 +232,7 @@ Constraint {
   { Name BlochY;
     Case {
       { Region SurfBlochYp; Type LinkCplx ; RegionRef SurfBlochYm;
-      Coefficient dephY[]; Function Vector[$X,$Y-period_y,$Z] ; }
+      Coefficient dephY[]; Function Vector[$X-dys,$Y-dyc,$Z] ; }
     }
   }
 }
@@ -330,6 +335,13 @@ PostProcessing {
   { Name postpro_helmholtz_vector; NameOfFormulation helmholtz_vector; NameOfSystem M;
     Quantity {
       { Name u      ; Value { Local { [ {u}       ]; In Omega; Jacobian JVol; } } }
+      { Name uper   ; Value { Local { [       {u} *Exp[-I[]*(skx1[]*X[]+sky1[]*Y[])]  ]; In Omega; Jacobian JVol; } } }
+      { Name uperx  ; Value { Local { [ CompX[{u} *Exp[-I[]*(skx1[]*X[]+sky1[]*Y[])]] ]; In Omega; Jacobian JVol; } } }
+      { Name upery  ; Value { Local { [ CompY[{u} *Exp[-I[]*(skx1[]*X[]+sky1[]*Y[])]] ]; In Omega; Jacobian JVol; } } }
+      { Name uperz  ; Value { Local { [ CompZ[{u} *Exp[-I[]*(skx1[]*X[]+sky1[]*Y[])]] ]; In Omega; Jacobian JVol; } } }
+      { 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 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; } } }
@@ -341,12 +353,15 @@ PostProcessing {
       { Name epsr_xx; Value { Local { [  CompXX[epsr[]] ]; In Omega; Jacobian JVol; } } }
       { Name Damp_pml_top; Value { Local { [Damp_pml_top[]  ]; In Omega; Jacobian JVol; } } }
       { Name Poy_tot; Value { Local { [ 0.5*Re[Cross[{u}+E1[] , Conj[ H1[]-I[]/(mur[]*mu0*om0)*{Curl u}]]] ]; In Omega; Jacobian JVol; } } }
+      { Name Poy_ref; Value { Local { [ 0.5*Re[Cross[{u}+E1d[], Conj[H1d[]-I[]/(mur[]*mu0*om0)*{Curl u}]]] ]; In Omega; Jacobian JVol; } } }
+      { Name Poy_inc; Value { Local { [ 0.5*Re[Cross[    Ei[] , Conj[Hi[]]]                              ] ]; In Omega; Jacobian JVol; } } }
       { Name lambda_step   ; Value { Local { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
 
       For k In {2:6}
-        { Name Abs_L~{k} ; Value { Integral { [ ep0*om0 * 0.5*Im[CompXX[epsr[]]]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*period_y) ] ; In L~{k} ; Integration I1 ; Jacobian JVol ; } } }
+        { Name Abs_L~{k} ; Value { Integral { [ ep0*om0 * 0.5*Im[CompXX[epsr[]]]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*dyc) ] ; In L~{k} ; Integration I1 ; Jacobian JVol ; } } }
       EndFor
-      { Name Abs_scat ; Value { Integral { [ ep0*om0 * 0.5*Im[CompXX[epsr[]]]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*period_y) ] ; In Scat ; Integration I1 ; Jacobian JVol ; } } }
+      { Name Abs_scat ; Value { Integral { [ ep0*om0 * 0.5*Im[CompXX[epsr[]]]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*dyc) ] ; In Scat ; Integration I1 ; Jacobian JVol ; } } }
+      { Name Abs_scat2; Value { Integral { [ ep0*om0 * 0.5*Im[CompXX[epsr[]]]*(SquNorm[{u}+E1[]]) ] ; In Scat ; Integration I1 ; Jacobian JVol ; } } }
 
       For i In {0:Nb_ordre-1}
         For j In {0:Nb_ordre-1}
@@ -355,11 +370,11 @@ PostProcessing {
           { 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 eff_t~{i}~{j}   ; Value { Term{Type Global; [
-                 1/(gammat~{i}~{j}[]*-CompZ[k1[]]) *((gammat~{i}~{j}[]^2+alpha~{i}[]^2)*SquNorm[$int_x_t~{i}~{j}]+
+                 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 ; } } }
           { Name eff_r~{i}~{j}   ; Value { Term{Type Global; [
-                 1/(gammar~{i}~{j}[]*-CompZ[k1[]])*((gammar~{i}~{j}[]^2+alpha~{i}[]^2)*SquNorm[$int_x_r~{i}~{j}]+
+                 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 ; } } }
           { Name numbering_ij~{i}~{j}   ; Value { Term{Type Global; [Vector[i-Nmax,j-Nmax,0]] ; In SurfIntBot ; } } }
@@ -423,6 +438,41 @@ PostOperation {
       // Print [ u , OnElementsOf Omega, File StrCat[myDir,"u.pos"]];
       // Print [ E1 , OnElementsOf Omega, File StrCat[myDir,"E1.pos"]];
       // Print [ Htotx , OnElementsOf Omega, File StrCat[myDir,"Htotx.pos"]];
+      // Print [ uper   , OnElementsOf Omega, File StrCat[myDir,"uper.pos"]];
+      // Print [ uperx  , OnElementsOf Omega, File StrCat[myDir,"uperx.pos"]];
+      // 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 [ 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} }
+                                  {npts_interpX-1,npts_interpY-1} , File StrCat[myDir,"Poy_tot.pos"]];
+      Print [ Poy_ref , OnPlane { {0.5*(-period_x-dys), -dyc/2, hh_L_1+thick_L_1/2}
+                                  {0.5*( period_x-dys), -dyc/2, hh_L_1+thick_L_1/2}
+                                  {0.5*(-period_x+dys),  dyc/2, hh_L_1+thick_L_1/2} }
+                                  {npts_interpX-1,npts_interpY-1} , File StrCat[myDir,"Poy_ref.pos"]];
+      Print [ Poy_inc , OnPlane { {0.5*(-period_x-dys), -dyc/2, hh_L_1+thick_L_1/2}
+                                  {0.5*( period_x-dys), -dyc/2, hh_L_1+thick_L_1/2}
+                                  {0.5*(-period_x+dys),  dyc/2, hh_L_1+thick_L_1/2} }
+                                  {npts_interpX-1,npts_interpY-1} , File StrCat[myDir,"Poy_inc.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} }
+                                  {npts_interpX-1,npts_interpY-1} , File StrCat[myDir,"Poy_tot_gd.pos"], Format Table];
+      Print [ Poy_ref , OnPlane { {0.5*(-period_x-dys), -dyc/2, hh_L_1+thick_L_1/2}
+                                  {0.5*( period_x-dys), -dyc/2, hh_L_1+thick_L_1/2}
+                                  {0.5*(-period_x+dys),  dyc/2, hh_L_1+thick_L_1/2} }
+                                  {npts_interpX-1,npts_interpY-1} , File StrCat[myDir,"Poy_ref_gd.pos"], Format Table];
+      Print [ Poy_inc , OnPlane { {0.5*(-period_x-dys), -dyc/2, hh_L_1+thick_L_1/2}
+                                  {0.5*( period_x-dys), -dyc/2, hh_L_1+thick_L_1/2}
+                                  {0.5*(-period_x+dys),  dyc/2, hh_L_1+thick_L_1/2} }
+                                  {npts_interpX-1,npts_interpY-1} , File StrCat[myDir,"Poy_inc_gd.pos"], Format Table];
+      Print[ Abs_scat2[Scat]  , OnGlobal, File > StrCat[myDir,"temp-Q_scat2.txt"], Format Table ];
+      //////////// END DEBUG  /////////
+
       For k In {2:6}
         Print[ Abs_L~{k}[L~{k}], OnGlobal, File > StrCat[myDir,Sprintf["temp-Q_L_%g.txt",k]], Format Table ];
       EndFor

DiffractionGratings/grating3D_postplot.py

 import numpy as np
 import matplotlib.pyplot as pl
 import sys
+def dtrap_poy(fname_in,nx,ny):
+    poy_data = np.loadtxt(fname_in)
+    x_export2D = poy_data[:,2].reshape((nx,ny))
+    y_export2D = poy_data[:,3].reshape((nx,ny))
+    poy_y_grid_re  = poy_data[:,10].reshape((nx,ny))
+    temp=np.trapz(poy_y_grid_re,x_export2D[:,0],axis=0)
+    return np.trapz(temp,y_export2D[0,:]) #[x_export2D,y_export2D,poy_y_grid_re] #
+
 myDir = sys.argv[1]
+
+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)
+Ascat2      = np.loadtxt(myDir+'/temp-Q_scat2.txt')[1]
+# poy_data = np.loadtxt('res3D/Poy_inc_gd.pos')
+# x_export2D = poy_data[:,2].reshape((50,50))
+# y_export2D = poy_data[:,3].reshape((50,50))
+# poy_y_grid_re  = poy_data[:,10].reshape((50,50))
+# temp=np.trapz(poy_y_grid_re,x_export2D[:,0],axis=0)
+# print(np.trapz(temp,y_export2D[0,:]))
+
 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]
 Q = [np.loadtxt(myDir+'/temp-Q_L_%g.txt'%k,ndmin=2)[:,1] for k in range(2,7)]
@@ -32,7 +52,11 @@ if myDir[6:]=='solarcell':
     pl.ylabel('fraction of incident energy')
     pl.savefig('solar_balance.pdf')
 else:
-    print('Rtot',Rnm.real.sum())
-    print('Ttot',Tnm.real.sum())
-    print('Atot',Q.sum())
-    print('TOT ',TOT)
+    print('Rtot ',Rnm.real.sum())
+    print('Rtot2',-intpoyz_ref/intpoyz_inc)
+    print('Ttot ',Tnm.real.sum())
+    print('Ttot2',intpoyz_tot/intpoyz_inc)
+    print('Atot ',Q.sum())
+    print('Atot2 ',Ascat2/-intpoyz_inc)
+    print('TOT   ',TOT)
+    print('TOT2  ',(-Ascat2+intpoyz_tot-intpoyz_ref)/intpoyz_inc)