From 5aa612e3be8aae5bf2959b3c9a49cd60a681d7da Mon Sep 17 00:00:00 2001
From: Guillaume Demesy <guillaume.demesy@fresnel.fr>
Date: Fri, 15 Nov 2019 15:20:10 +0100
Subject: [PATCH] better def of E in TM case
---
DiffractionGratings/grating2D.pro | 156 +++++++++---------
DiffractionGratings/grating2D_data.geo | 4 +-
.../grating2D_data_AnisotropicGrating.geo | 4 +-
.../grating2D_data_LamellarGrating.geo | 4 +-
.../grating2D_data_PhotonicCrystalSlab.geo | 4 +-
.../grating2D_data_ResonantGrating.geo | 4 +-
DiffractionGratings/grating2D_postplot.py | 6 +-
7 files changed, 89 insertions(+), 93 deletions(-)
diff --git a/DiffractionGratings/grating2D.pro b/DiffractionGratings/grating2D.pro
index 0dbed3f..f80a992 100644
--- a/DiffractionGratings/grating2D.pro
+++ b/DiffractionGratings/grating2D.pro
@@ -113,9 +113,8 @@ Function{
epsr1_re[] = epsr_re_dom_6[];
epsr1_im[] = epsr_im_dom_6[];
- Freq = cel/lambda0;
- omega0 = 2.*Pi*cel/lambda0;
- k0 = 2.*Pi/lambda0;
+ om0 = 2.*Pi*cel/lambda0;
+ k0 = 2.*Pi/lambda0;
epsr1[] = Complex[epsr1_re[],epsr1_im[]];
epsr2[] = Complex[epsr2_re[],epsr2_im[]];
n1[] = Sqrt[epsr1[]];
@@ -133,12 +132,12 @@ Function{
If (flag_Hparallel==1)
r[] = (CompY[k1[]]*epsr2[]-CompY[k2[]]*epsr1[])/(CompY[k1[]]*epsr2[]+CompY[k2[]]*epsr1[]);
t[] = (2.*CompY[k1[]]*epsr2[])/(CompY[k1[]]*epsr2[]+CompY[k2[]]*epsr1[]);
- Pinc[] = 0.5*A^2*Sqrt[mu0/(epsilon0*epsr1_re[])] * Cos[theta_deg*deg2rad];
+ Pinc[] = 0.5*A^2*Sqrt[mu0/(ep0*epsr1_re[])] * Cos[theta_deg*deg2rad];
EndIf
If (flag_Hparallel==0)
r[] = (CompY[k1[]]-CompY[k2[]])/(CompY[k1[]]+CompY[k2[]]);
t[] = 2.*CompY[k1[]] /(CompY[k1[]]+CompY[k2[]]);
- Pinc[] = 0.5*A^2*Sqrt[epsilon0*epsr1_re[]/mu0] * Cos[theta_deg*deg2rad];
+ Pinc[] = 0.5*A^2*Sqrt[ep0*epsr1_re[]/mu0] * Cos[theta_deg*deg2rad];
EndIf
BlochX_phase_shift[] = Exp[I[]*alpha[]*d];
@@ -228,19 +227,13 @@ Function{
u1d[] = ur[]+ut[];
If (flag_Hparallel==1)
- Exi[] = CompY[k1[] ]*ui[] / (omega0*epsilon0*CompXX[epsilonr_annex[]]);
- Exr[] = CompY[k1r[]]*ur[] / (omega0*epsilon0*CompXX[epsilonr_annex[]]);
- Ext[] = CompY[k2[] ]*ut[] / (omega0*epsilon0*CompXX[epsilonr_annex[]]);
- Eyi[] = -CompX[k1[] ]*ui[] / (omega0*epsilon0*CompXX[epsilonr_annex[]]);
- Eyr[] = -CompX[k1r[]]*ur[] / (omega0*epsilon0*CompXX[epsilonr_annex[]]);
- Eyt[] = -CompX[k2[] ]*ut[] / (omega0*epsilon0*CompXX[epsilonr_annex[]]);
- Ex1[Omega_top] = Exi[]+Exr[];
- Ey1[Omega_top] = Eyi[]+Eyr[];
- Ex1[Omega_bot] = Ext[];
- Ey1[Omega_bot] = Eyt[];
- Ex1[Omega_pml] = 0;
- Ey1[Omega_pml] = 0;
- E1[] = Vector[Ex1[],Ey1[],0];
+ E1i[] = -1/(om0*ep0*epsilonr_annex[]) * k1[] /\ Vector[0,0,ui[]];
+ E1r[] = -1/(om0*ep0*epsilonr_annex[]) * k1r[] /\ Vector[0,0,ur[]];
+ E1t[] = -1/(om0*ep0*epsilonr_annex[]) * k2[] /\ Vector[0,0,ut[]];
+ E1[] = E1i[]+E1r[]+E1t[];
+ Ex1[] = CompX[E1[]];
+ Ey1[] = CompY[E1[]];
+
detepst[] = CompXX[epsilonr[]]*CompYY[epsilonr[]]-CompXY[epsilonr[]]*Conj[CompYX[epsilonr[]]];
detepst_annex[] = CompXX[epsilonr_annex[]]*CompYY[epsilonr_annex[]]-CompXY[epsilonr_annex[]]*Conj[CompYX[epsilonr_annex[]]];
xsi[] = Transpose[epsilonr[]]/detepst[];
@@ -250,6 +243,7 @@ Function{
source_xsi_i[] = (xsi[]-xsi_annex[]) * -k1r[] * I[] * ur[];
source_chi_r[] = 0;
source_chi_i[] = 0;
+
Else
detmut[] = CompXX[mur[]]*CompYY[mur[]]-CompXY[mur[]]*Conj[CompYX[mur[]]];
xsi[] = Transpose[mur[]]/detmut[];
@@ -354,7 +348,7 @@ Formulation {
Resolution {
{ Name helmoltz_scalar;
System {
- { Name S; NameOfFormulation helmoltz_scalar; Type ComplexValue; Frequency Freq;}
+ { Name S; NameOfFormulation helmoltz_scalar; Type ComplexValue; }
}
Operation {
CreateDir[Str[myDir]];
@@ -376,6 +370,7 @@ PostProcessing {
{ Name Hz_diff ; Value { Local { [ {u2d}+u1d[] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_tot ; Value { Local { [ {u2d}+u1[] ]; In Omega; Jacobian JVol; } } }
{ Name NormHz_tot ; Value { Local { [ Norm[{u2d}+u1[]] ]; In Omega; Jacobian JVol; } } }
+ { Name E1 ; Value { Local { [ E1[] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_totp1 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 1*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_totp2 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 2*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_totp3 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 3*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
@@ -385,8 +380,6 @@ PostProcessing {
{ Name Hz_totm3 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-3*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_totm4 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-4*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name boundary ; Value { Local { [ bndCol[] ] ; In Plot_bnd ; Jacobian JVol ; } } }
-
- { Name E1 ; Value { Local { [ E1[] ]; In Omega; Jacobian JVol; } } }
For i In {0:2*nb_orders}
{ Name s_r~{i} ; Value { Integral{ [ expmialpha~{i}[] * ({u2d}+u1d[])/d ] ; In SurfCutSuper1 ; Jacobian JSur ; Integration Int_1 ; } } }
@@ -399,13 +392,13 @@ PostProcessing {
{ Name eff_t~{i} ; Value { Term{ Type Global; [ SquNorm[#(2*nb_orders+1+i)]*(beta_subs~{i}[]/beta1[])*(epsr1[]/epsr2[]) ] ; In SurfCutSubs1 ; } } }
EndFor
For i In {0:N_rods-1:1}
- { Name Q_rod~{i} ; Value { Integral { [ 0.5 * epsilon0*omega0*Fabs[epsr_rods_im[]] * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[ CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ey1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] ) / (Pinc[]*d) ] ; In rod~{i} ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Q_rod~{i} ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_rods_im[]] * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ey1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] ) / (Pinc[]*d) ] ; In rod~{i} ; Integration Int_1 ; Jacobian JVol ; } } }
EndFor
- { Name Q_subs ; Value { Integral { [ 0.5 * epsilon0*omega0*Fabs[epsr2_im[] ] * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In sub ; Integration Int_1 ; Jacobian JVol ; } } }
- { Name Q_rod_out ; Value { Integral { [ 0.5 * epsilon0*omega0*Fabs[epsr_rod_out_im[]] * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In rod_out ; Integration Int_1 ; Jacobian JVol ; } } }
- { Name Q_layer_dep ; Value { Integral { [ 0.5 * epsilon0*omega0*Fabs[epsr_layer_dep_im[]] * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In layer_dep ; Integration Int_1 ; Jacobian JVol ; } } }
- { Name Q_layer_cov ; Value { Integral { [ 0.5 * epsilon0*omega0*Fabs[epsr_layer_cov_im[]] * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In layer_cov ; Integration Int_1 ; Jacobian JVol ; } } }
- { Name Q_tot ; Value { Integral { [ 0.5 * epsilon0*omega0*Fabs[Im[CompZZ[epsilonr[]]]] * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In Plot_domain ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Q_subs ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr2_im[] ] * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In sub ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Q_rod_out ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_rod_out_im[]] * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In rod_out ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Q_layer_dep ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_layer_dep_im[]] * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In layer_dep ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Q_layer_cov ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_layer_cov_im[]] * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In layer_cov ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Q_tot ; Value { Integral { [ 0.5 * ep0*om0*Fabs[Im[CompZZ[epsilonr[]]]] * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In Plot_domain ; Integration Int_1 ; Jacobian JVol ; } } }
{ Name lambda_step ; Value { Local { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
EndIf
If (flag_Hparallel==0)
@@ -427,22 +420,22 @@ PostProcessing {
{ Name s_t~{i} ; Value { Integral{ [ expmialpha~{i}[] * ({u2d}+u1d[])/d ] ; In SurfCutSubs1 ; Jacobian JSur ; Integration Int_1 ; } } }
{ Name order_t_angle~{i} ; Value { Local{ [-Atan2[Re[alpha~{i}[]],Re[beta_subs~{i}[]]]/deg2rad ] ; In Omega; Jacobian JVol; } } }
{ Name order_r_angle~{i} ; Value { Local{ [ Atan2[Re[alpha~{i}[]],Re[beta_super~{i}[]]]/deg2rad ] ; In Omega; Jacobian JVol; } } }
- EndFor
+ EndFor
For i In {0:2*nb_orders}
{ Name eff_r~{i} ; Value { Term{ Type Global; [ SquNorm[#i]*beta_super~{i}[]/beta1[] ] ; In SurfCutSuper1 ; } } }
{ Name eff_t~{i} ; Value { Term{ Type Global; [ SquNorm[#(2*nb_orders+1+i)]*(beta_subs~{i}[]/beta1[])] ; In SurfCutSubs1 ; } } }
EndFor
For i In {0:N_rods-1:1}
- { Name Q_rod~{i} ; Value { Integral { [0.5 * epsilon0*omega0*Fabs[epsr_rods_im[]] * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In rod~{i} ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Q_rod~{i} ; Value { Integral { [0.5 * ep0*om0*Fabs[epsr_rods_im[]] * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In rod~{i} ; Integration Int_1 ; Jacobian JVol ; } } }
EndFor
- { Name Q_subs ; Value { Integral { [ 0.5 * epsilon0*omega0*Fabs[epsr2_im[]] *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In sub ; Integration Int_1 ; Jacobian JVol ; } } }
- { Name Q_rod_out ; Value { Integral { [ 0.5 * epsilon0*omega0*Fabs[epsr_rod_out_im[]] *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In rod_out ; Integration Int_1 ; Jacobian JVol ; } } }
- { Name Q_layer_dep ; Value { Integral { [ 0.5 * epsilon0*omega0*Fabs[epsr_layer_dep_im[]] *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In layer_dep ; Integration Int_1 ; Jacobian JVol ; } } }
- { Name Q_layer_cov ; Value { Integral { [ 0.5 * epsilon0*omega0*Fabs[epsr_layer_cov_im[]] *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In layer_cov ; Integration Int_1 ; Jacobian JVol ; } } }
- { Name Q_tot ; Value { Integral { [ 0.5 * epsilon0*omega0*Fabs[Im[CompXX[epsilonr[]]]]*(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In Plot_domain ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Q_subs ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr2_im[]] *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In sub ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Q_rod_out ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_rod_out_im[]] *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In rod_out ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Q_layer_dep ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_layer_dep_im[]] *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In layer_dep ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Q_layer_cov ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_layer_cov_im[]] *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In layer_cov ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Q_tot ; Value { Integral { [ 0.5 * ep0*om0*Fabs[Im[CompXX[epsilonr[]]]]*(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In Plot_domain ; Integration Int_1 ; Jacobian JVol ; } } }
{ Name lambda_step ; Value { Local { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
EndIf
- }
+ }
}
}
@@ -467,11 +460,12 @@ PostOperation {
For i In {0:N_rods-1:1}
Print[ Q_rod~{i}[rod~{i}] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, Sprintf("absorption-Q_rod_%g.txt", i+1) ]];
EndFor
- Print[ Q_subs[sub] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_subs.txt" ]];
- Print[ Q_rod_out[rod_out] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_rod_out.txt"]];
+ Print[ Q_subs[sub] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_subs.txt"]];
+ Print[ Q_rod_out[rod_out] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_rod_out.txt"]];
Print[ Q_layer_dep[layer_dep] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_dep.txt"]];
Print[ Q_layer_cov[layer_cov] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_cov.txt"]];
Print[ Hz_tot , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm)];
+ Print[ E1 , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("E1%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm)];
Echo[ Str["For i In {0:2}",
" View[i].LineWidth = 4;View[i].ColormapNumber = 15;View[0].AxesFormatX = '%.1f';",
"EndFor"], File StrCat[myDir,"tmp0.geo"]] ;
@@ -539,53 +533,55 @@ PostOperation {
" View[i].LineWidth = 4;View[i].ColormapNumber = 15;View[0].AxesFormatX = '%.1f';",
"EndFor"], File StrCat[myDir,"tmp0.geo"]] ;
If(multiplot)
- Echo[ Str["For i In {PostProcessing.NbViews-1:0:-1}",
- " If(!StrCmp(View[i].Name, 'boundary') || !StrCmp(View[i].Name, 'boundary_Combine'))",
- " Delete View[i];",
- " EndIf",
- "EndFor"], File StrCat[myDir,"tmp1.geo"]] ;
- Print [ Ez_totp1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_2.pos", lambda0/nm)], ChangeOfCoordinates {$X+1*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
- Print [ Ez_totm1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_3.pos", lambda0/nm)], ChangeOfCoordinates {$X-1*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
- Print [ Ez_totp2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_4.pos", lambda0/nm)], ChangeOfCoordinates {$X+2*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
- Print [ Ez_totm2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_5.pos", lambda0/nm)], ChangeOfCoordinates {$X-2*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
- Print [ Ez_totp3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_6.pos", lambda0/nm)], ChangeOfCoordinates {$X+3*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
- Print [ Ez_totm3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_7.pos", lambda0/nm)], ChangeOfCoordinates {$X-3*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
- Print [ Ez_totp4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_8.pos", lambda0/nm)], ChangeOfCoordinates {$X+4*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
- Print [ Ez_totm4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_9.pos", lambda0/nm)], ChangeOfCoordinates {$X-4*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
- Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo"]] ;
- Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary1.pos"], ChangeOfCoordinates {$X+0*d,$Y,$Z} ] ;
- Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary2.pos"], ChangeOfCoordinates {$X+1*d,$Y,$Z} ] ;
- Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary3.pos"], ChangeOfCoordinates {$X-1*d,$Y,$Z} ] ;
- Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary4.pos"], ChangeOfCoordinates {$X+2*d,$Y,$Z} ] ;
- Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary5.pos"], ChangeOfCoordinates {$X-2*d,$Y,$Z} ] ;
- Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary6.pos"], ChangeOfCoordinates {$X+3*d,$Y,$Z} ] ;
- Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary7.pos"], ChangeOfCoordinates {$X-3*d,$Y,$Z} ] ;
- Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary8.pos"], ChangeOfCoordinates {$X+4*d,$Y,$Z} ] ;
- Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary9.pos"], ChangeOfCoordinates {$X-4*d,$Y,$Z} ];
- Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo" ]] ;
- Echo[ Str["Hide {",
- "Point{1,2,7,8,9,10,20,22};",
- "Line{1,7,8,9,10,30,32,34,2,3,4,5,6,12,16,20,24,28};",
- "Surface{36,48};}",
- "Geometry.Color.Lines = {0,0,0};",
- "l=PostProcessing.NbViews-1; View[l].ColorTable={Black}; ",
- "View[l-1].Visible=1; View[l-1].ShowScale=0;",
- "View[l].ShowScale=0; View[l].LineWidth=1.5; View[l].LineType=1;Geometry.LineWidth=0;"],
- File StrCat[myDir,"tmp3.geo" ]] ;
- EndIf
+ Echo[ Str["For i In {PostProcessing.NbViews-1:0:-1}",
+ " If(!StrCmp(View[i].Name, 'boundary') || !StrCmp(View[i].Name, 'boundary_Combine'))",
+ " Delete View[i];",
+ " EndIf",
+ "EndFor"], File StrCat[myDir,"tmp1.geo"]];
+ Print [ Ez_totp1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_2.pos", lambda0/nm)], ChangeOfCoordinates {$X+1*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+ Print [ Ez_totm1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_3.pos", lambda0/nm)], ChangeOfCoordinates {$X-1*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+ Print [ Ez_totp2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_4.pos", lambda0/nm)], ChangeOfCoordinates {$X+2*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+ Print [ Ez_totm2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_5.pos", lambda0/nm)], ChangeOfCoordinates {$X-2*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+ Print [ Ez_totp3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_6.pos", lambda0/nm)], ChangeOfCoordinates {$X+3*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+ Print [ Ez_totm3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_7.pos", lambda0/nm)], ChangeOfCoordinates {$X-3*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+ Print [ Ez_totp4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_8.pos", lambda0/nm)], ChangeOfCoordinates {$X+4*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+ Print [ Ez_totm4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_9.pos", lambda0/nm)], ChangeOfCoordinates {$X-4*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+ Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo"]] ;
+ Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary1.pos"], ChangeOfCoordinates {$X+0*d,$Y,$Z} ] ;
+ Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary2.pos"], ChangeOfCoordinates {$X+1*d,$Y,$Z} ] ;
+ Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary3.pos"], ChangeOfCoordinates {$X-1*d,$Y,$Z} ] ;
+ Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary4.pos"], ChangeOfCoordinates {$X+2*d,$Y,$Z} ] ;
+ Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary5.pos"], ChangeOfCoordinates {$X-2*d,$Y,$Z} ] ;
+ Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary6.pos"], ChangeOfCoordinates {$X+3*d,$Y,$Z} ] ;
+ Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary7.pos"], ChangeOfCoordinates {$X-3*d,$Y,$Z} ] ;
+ Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary8.pos"], ChangeOfCoordinates {$X+4*d,$Y,$Z} ] ;
+ Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary9.pos"], ChangeOfCoordinates {$X-4*d,$Y,$Z} ];
+ Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo" ]] ;
+ Echo[ Str["Hide {",
+ "Point{1,2,7,8,9,10,20,22};",
+ "Line{1,7,8,9,10,30,32,34,2,3,4,5,6,12,16,20,24,28};",
+ "Surface{36,48};}",
+ "Geometry.Color.Lines = {0,0,0};",
+ "l=PostProcessing.NbViews-1; View[l].ColorTable={Black}; ",
+ "View[l-1].Visible=1; View[l-1].ShowScale=0;",
+ "View[l].ShowScale=0; View[l].LineWidth=1.5; View[l].LineType=1;Geometry.LineWidth=0;"],
+ File StrCat[myDir,"tmp3.geo" ]];
+ EndIf
EndIf
}
}
}
+
DefineConstant[
- R_ = {"helmoltz_scalar", Name "GetDP/1ResolutionChoices", Visible 1},
- C_ = {"-solve -pos -petsc_prealloc 100 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
- P_ = {"postop_energy", Name "GetDP/2PostOperationChoices", Visible 1}];
-
+ R_ = {"helmoltz_scalar", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {"-solve -pos -petsc_prealloc 100 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {"postop_energy", Name "GetDP/2PostOperationChoices", Visible 1}
+];
+
If(plotRTgraphs)
-DefineConstant[
- refl_ = {0, Name "GetDP/R0", ReadOnly 1, Graph "02000000", Visible 1},
- abs_ = {0, Name "GetDP/total absorption", ReadOnly 1, Graph "00000002", Visible 1},
- trans_ = {0, Name "GetDP/T0", ReadOnly 1, Graph "000000000002", Visible 1}
- ];
+ DefineConstant[
+ refl_ = {0, Name "GetDP/R0", ReadOnly 1, Graph "02000000", Visible 1},
+ abs_ = {0, Name "GetDP/total absorption", ReadOnly 1, Graph "00000002", Visible 1},
+ trans_ = {0, Name "GetDP/T0", ReadOnly 1, Graph "000000000002", Visible 1}
+ ];
EndIf
diff --git a/DiffractionGratings/grating2D_data.geo b/DiffractionGratings/grating2D_data.geo
index 0a11dac..61f5259 100644
--- a/DiffractionGratings/grating2D_data.geo
+++ b/DiffractionGratings/grating2D_data.geo
@@ -4,9 +4,9 @@
///////////////////////////////////
nm = 1000.;
-epsilon0 = 8.854187817e-3*nm;
+ep0 = 8.854187817e-3*nm;
mu0 = 400.*Pi*nm;
-cel = 1.0/(Sqrt[epsilon0 * mu0]);
+cel = 1.0/(Sqrt[ep0 * mu0]);
deg2rad = Pi/180;
pp0 = "1Geometry/0";
diff --git a/DiffractionGratings/grating2D_data_AnisotropicGrating.geo b/DiffractionGratings/grating2D_data_AnisotropicGrating.geo
index 5296621..0036f31 100644
--- a/DiffractionGratings/grating2D_data_AnisotropicGrating.geo
+++ b/DiffractionGratings/grating2D_data_AnisotropicGrating.geo
@@ -4,9 +4,9 @@
////////////////////////////////////////////
nm = 1000.;
-epsilon0 = 8.854187817e-3*nm;
+ep0 = 8.854187817e-3*nm;
mu0 = 400.*Pi*nm;
-cel = 1.0/(Sqrt[epsilon0 * mu0]);
+cel = 1.0/(Sqrt[ep0 * mu0]);
deg2rad = Pi/180;
diff --git a/DiffractionGratings/grating2D_data_LamellarGrating.geo b/DiffractionGratings/grating2D_data_LamellarGrating.geo
index d4436bf..d31b0d8 100644
--- a/DiffractionGratings/grating2D_data_LamellarGrating.geo
+++ b/DiffractionGratings/grating2D_data_LamellarGrating.geo
@@ -4,9 +4,9 @@
/////////////////////////////////////////
nm = 1000.;
-epsilon0 = 8.854187817e-3*nm;
+ep0 = 8.854187817e-3*nm;
mu0 = 400.*Pi*nm;
-cel = 1.0/(Sqrt[epsilon0 * mu0]);
+cel = 1.0/(Sqrt[ep0 * mu0]);
deg2rad = Pi/180;
diff --git a/DiffractionGratings/grating2D_data_PhotonicCrystalSlab.geo b/DiffractionGratings/grating2D_data_PhotonicCrystalSlab.geo
index 509c764..9618055 100644
--- a/DiffractionGratings/grating2D_data_PhotonicCrystalSlab.geo
+++ b/DiffractionGratings/grating2D_data_PhotonicCrystalSlab.geo
@@ -4,9 +4,9 @@
/////////////////////////////////////////////
nm = 1000.;
-epsilon0 = 8.854187817e-3*nm;
+ep0 = 8.854187817e-3*nm;
mu0 = 400.*Pi*nm;
-cel = 1.0/(Sqrt[epsilon0 * mu0]);
+cel = 1.0/(Sqrt[ep0 * mu0]);
deg2rad = Pi/180;
pp0 = "1Geometry/0";
diff --git a/DiffractionGratings/grating2D_data_ResonantGrating.geo b/DiffractionGratings/grating2D_data_ResonantGrating.geo
index 786c194..ce1b512 100644
--- a/DiffractionGratings/grating2D_data_ResonantGrating.geo
+++ b/DiffractionGratings/grating2D_data_ResonantGrating.geo
@@ -4,9 +4,9 @@
/////////////////////////////////////////
nm = 1000.;
-epsilon0 = 8.854187817e-3*nm;
+ep0 = 8.854187817e-3*nm;
mu0 = 400.*Pi*nm;
-cel = 1.0/(Sqrt[epsilon0 * mu0]);
+cel = 1.0/(Sqrt[ep0 * mu0]);
deg2rad = Pi/180;
pp0 = "1Geometry/0";
pp01 = "1Geometry/01Stack thicknesses/";
diff --git a/DiffractionGratings/grating2D_postplot.py b/DiffractionGratings/grating2D_postplot.py
index 9da09b3..8504011 100644
--- a/DiffractionGratings/grating2D_postplot.py
+++ b/DiffractionGratings/grating2D_postplot.py
@@ -14,9 +14,9 @@ pl.rc('legend',fontsize=20)
# pl.rc('figure', **{'autolayout': True})
#### scaling the problem (*10^12)
nm = 1.e3
-epsilon0 = 8.854187817e-3*nm
+ep0 = 8.854187817e-3*nm
mu0 = 400.*pi*nm
-cel = 1.0/(np.sqrt(epsilon0 * mu0))
+cel = 1.0/(np.sqrt(ep0 * mu0))
####
nb_orders = int(int(subprocess.check_output("ls ./run_results/efficiency_r_* | grep -c efficiency_r_", shell=True))/2)
nb_rods = int(int(subprocess.check_output("ls ./run_results/absorption-Q_rod_* | grep -c absorption-Q_rod_", shell=True))-1)
@@ -42,7 +42,7 @@ if len(np.loadtxt('./run_results/efficiency_r_0.txt').shape)==2:
A_rod_out = np.loadtxt('./run_results/absorption-Q_rod_out.txt')[:,1]
A_layer_cov = np.loadtxt('./run_results/absorption-Q_layer_cov.txt')[:,1]
A_layer_dep = np.loadtxt('./run_results/absorption-Q_layer_dep.txt')[:,1]
- A_sub = np.loadtxt('./run_results/absorption-Q_sub.txt')[:,1]
+ A_sub = np.loadtxt('./run_results/absorption-Q_subs.txt')[:,1]
pl.savez('last_run_RTA.npz',R0=R0,T0=T0,A=A)
--
GitLab