diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index 8ef65f80b4a2deb5f2e9d399ee80b187e8c4f852..fc14c85f99cf9f47410be02b6c8b4c68e5cc3aed 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/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
diff --git a/DiffractionGratings/grating3D_postplot.py b/DiffractionGratings/grating3D_postplot.py
index 2c3b472e801acd7984fb1cd15b492a6a92715755..ad1cedd326d66803526e4433e840b273f3b817cd 100644
--- a/DiffractionGratings/grating3D_postplot.py
+++ b/DiffractionGratings/grating3D_postplot.py
@@ -1,7 +1,27 @@
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)