From 18b5e01e5b5a5d296f8eb54d247ff58771620cf8 Mon Sep 17 00:00:00 2001
From: Guillaume Demesy <guillaume.demesy@fresnel.fr>
Date: Thu, 22 Sep 2022 22:26:48 +0200
Subject: [PATCH] polarimetric quantities for M-matrix computations
---
DiffractionGratings/grating3D.geo | 5 +
DiffractionGratings/grating3D.pro | 130 +++++++++++-------
.../grating3D_data_halfellipsoid.geo | 10 +-
DiffractionGratings/grating3D_postplot.py | 63 ++++++---
DiffractionGratings/grating3D_runall.sh | 7 +-
5 files changed, 141 insertions(+), 74 deletions(-)
diff --git a/DiffractionGratings/grating3D.geo b/DiffractionGratings/grating3D.geo
index 31934f7..69c2217 100644
--- a/DiffractionGratings/grating3D.geo
+++ b/DiffractionGratings/grating3D.geo
@@ -289,3 +289,8 @@ EndIf
// Mesh.SurfaceEdges = 0;
Mesh.VolumeEdges = 0;
Mesh.ElementOrder = og;
+If (og==2)
+ Mesh.HighOrderOptimize = 1;
+EndIf
+
+Solver.SocketName=Sprintf("socktest%g",socketid);
\ No newline at end of file
diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index 1a4741b..07461d4 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -53,14 +53,16 @@ Group {
L_1 = Region[{L_1_temp,SurfIntTop}];
L_6 = Region[{L_6_temp,SurfIntBot}];
- Omega = Region[{PMLbot,L_6,L_5,L_4,L_3,L_2,L_1,PMLtop,Scat}];
- Omega_nosource = Region[{PMLbot,L_6,L_1,PMLtop}];
- Omega_source = Region[{Scat,L_2,L_3,L_4,L_5}];
- Omega_super = Region[{Scat,L_1,L_2,L_3,L_4,L_5,PMLtop}];
- Omega_subs = Region[{L_6,PMLbot}];
- Omega_plot = Region[{L_6,L_5,L_4,L_3,L_2,L_1,Scat}];
-
- // SurfNeumann = Region[{SurfBlochXm,SurfBlochXp,SurfBlochYm,SurfBlochYp}];
+ Omega = Region[{PMLbot,L_6,L_5,L_4,L_3,L_2,L_1,PMLtop,Scat}];
+ Omega_nosource = Region[{PMLbot,L_6,L_1,PMLtop}];
+ Omega_source = Region[{Scat,L_2,L_3,L_4,L_5}];
+ Omega_super = Region[{Scat,L_1,L_2,L_3,L_4,L_5,PMLtop}];
+ Omega_subs = Region[{L_6,PMLbot}];
+ Omega_plot = Region[{L_6,L_5,L_4,L_3,L_2,L_1,Scat}];
+
+ // get normal comp of E field on integ surfaces
+ Gama = Region[{SurfIntBot,SurfIntTop}];
+ Tr = ElementsOf[Omega_plot, ConnectedTo Gama];
}
@@ -69,6 +71,9 @@ Function{
I[] = Complex[0,1];
zhat[] = Vector[0,0,1];
+ ispecular = Nmax;
+ jspecular = Nmax;
+
small_delta = 0.0*nm;
mu0 = 4*Pi*1.e2*nm;
ep0 = 8.854187817e-3*nm;
@@ -166,7 +171,9 @@ Function{
rp[] = (k1z[]*epsr2[]-k2z[]*epsr1[])/(k1z[]*epsr2[]+k2z[]*epsr1[]);
tp[] = (2.*k1z[]*epsr2[])/(k1z[]*epsr2[]+k2z[]*epsr1[]);
- spol[] = Vector[Sin[phi0],-Cos[phi0],0];
+ spol[] = Vector[Sin[phi0],-Cos[phi0],0];
+ ppol[] = (k1r[]/Norm[k1r[]]) /\ spol[];
+
AmplEis[] = spol[];
AmplErs[] = rs[]*spol[];
AmplEts[] = ts[]*spol[];
@@ -249,21 +256,9 @@ Constraint {
}
Jacobian {
- { Name JVol ;
- Case {
- { Region All ; Jacobian Vol ; }
- }
- }
- { Name JSur ;
- Case {
- { Region All ; Jacobian Sur ; }
- }
- }
- { Name JLin ;
- Case {
- { Region All ; Jacobian Lin ; }
- }
- }
+ { Name JVol ; Case {{ Region All ; Jacobian Vol ; }}}
+ { Name JSur ; Case {{ Region All ; Jacobian Sur ; }}}
+ { Name JLin ; Case {{ Region All ; Jacobian Lin ; }}}
}
Integration {
@@ -314,12 +309,19 @@ FunctionSpace {
{ NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
}
}
+ { Name L2_lambda; Type Form0;
+ BasisFunction{
+ { Name ln ; NameOfCoef lambda_n ; Function BF_Node ; Support Region[Gama]; Entity NodesOf[All];}
+ { Name ln2; NameOfCoef lambda_n2; Function BF_Node_2E; Support Region[Gama]; Entity EdgesOf[All];}
+ }
+ }
}
Formulation {
{ Name helmholtz_vector; Type FemEquation;
Quantity {
- { Name u; Type Local; NameOfSpace Hcurl;}
+ { Name u ; Type Local; NameOfSpace Hcurl; }
+ { Name uz ; Type Local; NameOfSpace L2_lambda; }
}
Equation{
Galerkin { [-1/mur[] * Dof{Curl u} , {Curl u}]; In Omega ; Jacobian JVol; Integration I1; }
@@ -329,6 +331,8 @@ Formulation {
Else
Galerkin { [source_surf_tot[] , {u}]; In SurfIntTop; Jacobian JVol; Integration I1; }
EndIf
+ Galerkin{ [ Dof{uz} , {uz}]; In Gama; Jacobian JSur; Integration I1;}
+ Galerkin{ [ Trace[Dof{u}, Tr]*Vector[0,0,-1] , {uz}]; In Gama; Jacobian JSur; Integration I1;}
}
}
}
@@ -350,6 +354,7 @@ PostProcessing {
{ Name postpro_helmholtz_vector; NameOfFormulation helmholtz_vector; NameOfSystem M;
Quantity {
{ Name u ; Value { Local { [ {u} ]; In Omega; Jacobian JVol; } } }
+ { Name uz ; Value { Local { [ {uz} ]; In Gama; Jacobian JSur; } } }
{ Name Damp_pml_top; Value { Local { [Damp_pml_top[] ]; In Omega; Jacobian JVol; } } }
{ Name epsr_xx; Value { Local { [ CompXX[epsr[]] ]; In Omega; Jacobian JVol; } } }
{ Name Poy_inc; Value { Local { [ 0.5*Re[Cross[ Ei[] , Conj[Hi[]]]] ]; In Omega; Jacobian JVol; } } }
@@ -379,34 +384,63 @@ PostProcessing {
For i In {0:Nb_ordre-1}
For j In {0:Nb_ordre-1}
If (FLAG_TOTAL==0)
- { 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 ; } } }
- Else
- { Name int_x_t~{i}~{j} ; Value { Integral { [ CompX[{u} ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
- { Name int_y_t~{i}~{j} ; Value { Integral { [ CompY[{u} ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
- { Name int_x_r~{i}~{j} ; Value { Integral { [ CompX[{u}-Ei[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
- { Name int_y_r~{i}~{j} ; Value { Integral { [ CompY[{u}-Ei[]]*expialphaxy~{i}~{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_z_t~{i}~{j} ; Value { Integral { [ ({uz}+CompZ[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 int_z_r~{i}~{j} ; Value { Integral { [({uz}+CompZ[E1d[]])*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
+ Else
+ { Name int_x_t~{i}~{j} ; Value { Integral { [ CompX[{u} ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
+ { Name int_y_t~{i}~{j} ; Value { Integral { [ CompY[{u} ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
+ { Name int_z_t~{i}~{j} ; Value { Integral { [ {uz} *expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
+ { Name int_x_r~{i}~{j} ; Value { Integral { [ CompX[{u}-Ei[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
+ { Name int_y_r~{i}~{j} ; Value { Integral { [ CompY[{u}-Ei[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
+ { Name int_z_r~{i}~{j} ; Value { Integral { [({uz}-CompZ[Ei[]])*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
EndIf
- { Name eff_t~{i}~{j} ; Value { Term{Type Global; [
+ { Name eff_t1~{i}~{j} ; Value { Term{ Type Global; [
1/(gammat~{i}~{j}[]*-k1z[]*Cos[xsi]^2) * ((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; [
+ { Name eff_r1~{i}~{j} ; Value { Term{ Type Global; [
1/(gammar~{i}~{j}[]*-k1z[]*Cos[xsi]^2) * ((gammar~{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 ; } } }
+ { Name eff_t2~{i}~{j} ; Value { Term{ Type Global; [
+ 1/(gammat~{i}~{j}[]*-k1z[]*Cos[xsi]^2) * (gammat~{i}~{j}[]^2 * SquNorm[$int_x_t~{i}~{j}]+
+ gammat~{i}~{j}[]^2 * SquNorm[$int_y_t~{i}~{j}]+
+ gammat~{i}~{j}[]^2 * SquNorm[$int_z_t~{i}~{j}] ) ] ; In SurfIntBot ; } } }
+ { Name eff_r2~{i}~{j} ; Value { Term{ Type Global; [
+ 1/(gammar~{i}~{j}[]*-k1z[]*Cos[xsi]^2) * (gammar~{i}~{j}[]^2 * SquNorm[$int_x_r~{i}~{j}]+
+ gammar~{i}~{j}[]^2 * SquNorm[$int_y_r~{i}~{j}]+
+ gammar~{i}~{j}[]^2 * SquNorm[$int_z_r~{i}~{j}] ) ] ; In SurfIntTop ; } } }
+ { Name numbering_ij~{i}~{j} ; Value { Term{ Type Global; [Vector[i-Nmax,j-Nmax,0]] ; In SurfIntBot ; } } }
EndFor
EndFor
- }
+ // Mmatrix computation : Retrieve the complex vector amplitude of the plane wave corresponding to the reflected specular order
+ // it is phase shifted by Exp[I[]*k1z[]*(thick_L_1+thick_L_2+thick_L_3)] because we measure it on SurfIntTop
+ // but it is better to have it at the surface on which the scatterer is relying (so that if there is no scatterer,
+ // it just corresponds to the usual definition of rs/rp for a simple diopter).
+ { Name er_specular ; Value { Term{ Type Global; [
+ Exp[I[]*k1z[]*(thick_L_1+thick_L_2+thick_L_3)] *
+ Vector[$int_x_r~{ispecular}~{jspecular},
+ $int_y_r~{ispecular}~{jspecular},
+ $int_z_r~{ispecular}~{jspecular}] ] ; In SurfIntTop ; } } }
+ // Project er_specular on the (s,p) basis
+ { Name rp ; Value { Term{ Type Global; [ ppol[] * $er_specular] ; In SurfIntTop ; } } }
+ { Name rs ; Value { Term{ Type Global; [ spol[] * $er_specular] ; In SurfIntTop ; } } }
+}
}
}
PostOperation {
{ Name postop_helmholtz_vector; NameOfPostProcessing postpro_helmholtz_vector ;
Operation {
+ Print [ uz , OnElementsOf Gama, File StrCat[myDir,"uz.pos"]];
+ Print [ u , OnPlane { {0.5*(-period_x-dys), -dyc/2,hh_L_6}
+ {0.5*( period_x-dys), -dyc/2,hh_L_6}
+ {0.5*(-period_x+dys), dyc/2,hh_L_6} }
+ {npts_checkpoyX-1,npts_checkpoyY-1} , File StrCat[myDir,"uz_check.pos"]];
If (FlagOutEscaFull==1)
If (Flag_interp_cubic==1)
Print [ u , OnBox { {-period_x/2,-period_y/2,hh_L_6-PML_bot} {period_x/2,-period_y/2,hh_L_6-PML_bot} {-period_x/2,period_y/2,hh_L_6-PML_bot} {-period_x/2,-period_y/2,hh_L_1+thick_L_1+PML_top} } {npts_interpX,npts_interpY,npts_interpZSca} , File StrCat[myDir,"u_grid.pos"], Name "u_grid"];
@@ -467,17 +501,24 @@ PostOperation {
For j In {0:Nb_ordre-1}
Print[ int_x_t~{i}~{j}[SurfIntBot], OnGlobal, StoreInVariable $int_x_t~{i}~{j}, Format Table];
Print[ int_y_t~{i}~{j}[SurfIntBot], OnGlobal, StoreInVariable $int_y_t~{i}~{j}, Format Table];
+ Print[ int_z_t~{i}~{j}[SurfIntBot], OnGlobal, StoreInVariable $int_z_t~{i}~{j}, Format Table];
Print[ int_x_r~{i}~{j}[SurfIntTop], OnGlobal, StoreInVariable $int_x_r~{i}~{j}, Format Table];
Print[ int_y_r~{i}~{j}[SurfIntTop], OnGlobal, StoreInVariable $int_y_r~{i}~{j}, Format Table];
+ Print[ int_z_r~{i}~{j}[SurfIntTop], OnGlobal, StoreInVariable $int_z_r~{i}~{j}, Format Table];
EndFor
EndFor
For i In {0:Nb_ordre-1}
For j In {0:Nb_ordre-1}
- Print[ eff_t~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir, "eff_t.txt"], Format Table ];
- Print[ eff_r~{i}~{j}[SurfIntTop], OnRegion SurfIntTop, File > StrCat[myDir, "eff_r.txt"], Format Table ];
+ Print[ eff_t1~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir, "eff_t1.txt"], Format Table ];
+ Print[ eff_r1~{i}~{j}[SurfIntTop], OnRegion SurfIntTop, File > StrCat[myDir, "eff_r1.txt"], Format Table ];
+ Print[ eff_t2~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir, "eff_t2.txt"], Format Table ];
+ Print[ eff_r2~{i}~{j}[SurfIntTop], OnRegion SurfIntTop, File > StrCat[myDir, "eff_r2.txt"], Format Table ];
EndFor
EndFor
+ Print[ er_specular[SurfIntTop] , OnRegion SurfIntTop, StoreInVariable $er_specular, Format Table];
+ Print[ rp[SurfIntTop], OnRegion SurfIntTop, File > StrCat[myDir,"rp.txt"], Format Table ];
+ Print[ rs[SurfIntTop], OnRegion SurfIntTop, File > StrCat[myDir,"rs.txt"], Format Table ];
For i In {0:Nb_ordre-1}
For j In {0:Nb_ordre-1}
Print[ numbering_ij~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir,"numbering_ij.txt"], Format Table ];
@@ -493,12 +534,3 @@ DefineConstant[
C_ = {"-solve -pos -petsc_prealloc 500 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps -ksp_error_if_not_converged", Name "GetDP/9ComputeCommand", Visible 1},
P_ = {"postop_helmholtz_vector", Name "GetDP/2PostOperationChoices", Visible 1}
];
-
-// Notes on possible mumps craches (NaN's when large nDOF)
-// -ksp_view
-// -ksp_error_if_not_converged
-// -ksp_monitor_true_residual
-// related to? https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2018-October/036377.html
-
-// N.B. pastix works with petsc 3.12
-// -pc_factor_mat_solver_type pastix
diff --git a/DiffractionGratings/grating3D_data_halfellipsoid.geo b/DiffractionGratings/grating3D_data_halfellipsoid.geo
index 9075c04..5095cab 100644
--- a/DiffractionGratings/grating3D_data_halfellipsoid.geo
+++ b/DiffractionGratings/grating3D_data_halfellipsoid.geo
@@ -6,8 +6,9 @@ pp4 = "4Layer Materials";
pp5 = "5Computational Parameters";
pp6 = "6Output";
DefineConstant[
- FLAG_TOTAL = {0 , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
- lambda0 = {495 , Name StrCat[pp1,"/1lambda0 [nm]"]},
+ FLAG_TOTAL = {0 , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
+ // lambda0 = {495 , Name StrCat[pp1,"/1lambda0 [nm]"]},
+ lambda0 = {400 , Min 400, Max 800, Step 200, Name StrCat[pp1, "0wavelength [nm]"] , Loop 1},
thetadeg = {40 , Name StrCat[pp1,"/2theta0 [deg]"]},
phideg = {36 , Name StrCat[pp1,"/3phi0 [deg]"]},
psideg = {72 , Name StrCat[pp1,"/4psi0 [deg]"]},
@@ -48,7 +49,7 @@ DefineConstant[
eps_re_L_6 = {4 , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
eps_im_L_6 = {0 , Name StrCat[pp4,"/layer 6: imag part of relative permittivity"]},
- og = {1 , Name StrCat[pp5,"/0geometrical order [-]"] , Choices {0="1",1="2"} },
+ og = {0 , Name StrCat[pp5,"/0geometrical order [-]"] , Choices {0="1",1="2"} },
oi = {1 , Name StrCat[pp5,"/0interpolation order [-]"], Choices {0="1",1="2"} },
paramaille = {8 , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
lc_scat = {10 , Name StrCat[pp5,"/2metal mesh size [nm]"]},
@@ -71,7 +72,8 @@ DefineConstant[
FlagOutPoyCut = { 1 , Name StrCat[pp6,"/5Output Poynting cuts?"] , Choices {0,1} },
FlagOutEtotFull = { 0 , Name StrCat[pp6,"/6Total Electric Field Full Output?"] , Choices {0,1} },
FlagOutEscaFull = { 0 , Name StrCat[pp6,"/7Scattered Electric Field Full Output?"] , Choices {0,1} },
- FlagOutPoyFull = { 0 , Name StrCat[pp6,"/8Poynting Full Output?"] , Choices {0,1} }
+ FlagOutPoyFull = { 0 , Name StrCat[pp6,"/8Poynting Full Output?"] , Choices {0,1} },
+ socketid = { 0 , Name StrCat[pp6,"/8socket ID"] }
];
lambda0 = nm * lambda0;
diff --git a/DiffractionGratings/grating3D_postplot.py b/DiffractionGratings/grating3D_postplot.py
index 80905d6..7fa2864 100644
--- a/DiffractionGratings/grating3D_postplot.py
+++ b/DiffractionGratings/grating3D_postplot.py
@@ -1,6 +1,23 @@
import numpy as np
import matplotlib.pyplot as pl
import sys
+
+termbg_black = '\033[107m'
+termbg_red = '\033[41m'
+termbg_green = '\033[42m'
+termbg_orange = '\033[43m'
+termbg_blue = '\033[44m'
+termbg_purple = '\033[45m'
+termbg_cyan = '\033[46m'
+termbg_lightgrey = '\033[47m'
+
+def colored_i(text): return "\033[38;2;{};{};{}m{} \033[38;2;255;255;255m".format(100 , 100, 100, text)
+def colored_R(text): return "\033[38;2;{};{};{}m{} \033[38;2;255;255;255m".format( 41, 84, 255, text)
+def colored_T(text): return "\033[38;2;{};{};{}m{} \033[38;2;255;255;255m".format(77 , 232, 159, text)
+def colored_A(text): return "\033[38;2;{};{};{}m{} \033[38;2;255;255;255m".format(250, 124, 112, text)
+def colored_TOT(text): return "\033[1m\033[38;2;{};{};{}m{} \033[38;2;255;255;255m\033[0m".format(245, 10, 10 , text)
+
+
def dtrap_poy(fname_in,nx,ny):
poy_data = np.loadtxt(fname_in)
x_export2D = poy_data[:,2].reshape((nx,ny))
@@ -14,23 +31,28 @@ myDir = sys.argv[1]
intpoyz_tot = abs(dtrap_poy(myDir+'/Poy_tot_gd.pos',50,50))
intpoyz_ref = abs(dtrap_poy(myDir+'/Poy_ref_gd.pos',50,50))
intpoyz_inc = abs(dtrap_poy(myDir+'/Poy_inc_gd.pos',50,50))
+R_poy = intpoyz_ref/intpoyz_inc
+T_poy = intpoyz_tot/intpoyz_inc
-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]
+R1nm = np.loadtxt(myDir+'/eff_r1.txt',ndmin=2)[:,1] + 1j*np.loadtxt(myDir+'/eff_r1.txt',ndmin=2)[:,2]
+T1nm = np.loadtxt(myDir+'/eff_t1.txt',ndmin=2)[:,1] + 1j*np.loadtxt(myDir+'/eff_t1.txt',ndmin=2)[:,2]
+R2nm = np.loadtxt(myDir+'/eff_r2.txt',ndmin=2)[:,1] + 1j*np.loadtxt(myDir+'/eff_r2.txt',ndmin=2)[:,2]
+T2nm = np.loadtxt(myDir+'/eff_t2.txt',ndmin=2)[:,1] + 1j*np.loadtxt(myDir+'/eff_t2.txt',ndmin=2)[:,2]
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)
-TOT = Rnm.real.sum()+Tnm.real.sum()+Q.sum()
+TOT1 = R1nm.real.sum()+T1nm.real.sum()+Q.sum()
+TOT2 = R2nm.real.sum()+T2nm.real.sum()+Q.sum()
if myDir[6:]=='solarcell':
print('cf pdf')
Nmax=2
tab_lambdas=np.loadtxt(myDir+'/temp_lambda_step.txt',ndmin=2)[:,8]
nb_lambdas=len(tab_lambdas)
- Rnm = Rnm.reshape((nb_lambdas,2*Nmax+1,2*Nmax+1))
- Tnm = Tnm.reshape((nb_lambdas,2*Nmax+1,2*Nmax+1))
- Rtot = [Rnm[i].real.sum() for i in range(nb_lambdas)]
- Ttot = [Tnm[i].real.sum() for i in range(nb_lambdas)]
+ R1nm = R1nm.reshape((nb_lambdas,2*Nmax+1,2*Nmax+1))
+ T1nm = T1nm.reshape((nb_lambdas,2*Nmax+1,2*Nmax+1))
+ Rtot = [R1nm[i].real.sum() for i in range(nb_lambdas)]
+ Ttot = [T1nm[i].real.sum() for i in range(nb_lambdas)]
Abs_rods = Q[-1]
Abs_ITO = Q[0]
Abs_subs = Q[2]+Q[3]+Q[4]+Ttot
@@ -54,14 +76,19 @@ elif myDir[6:]=='conv':
pl.legend()
pl.savefig('fig_convergence_absorption.pdf')
else:
- print('===> Computed from diffraction efficiencies')
- print('Rtot ',Rnm.real.sum())
- print('Ttot ',Tnm.real.sum())
- print('Atot ',Q.sum())
- print('TOT ',TOT)
- print('===> Computed from trapz on Poynting vector')
- print('(expected to be less precise)')
- print('Rtot2',intpoyz_ref/intpoyz_inc)
- print('Ttot2',intpoyz_tot/intpoyz_inc)
- print('Atot ',Q.sum())
- print('TOT2 ',Q.sum()+(intpoyz_tot+intpoyz_ref)/intpoyz_inc)
+ print(colored_i('===> Computed from diffraction efficiencies with tangential components only'))
+ print(colored_R('Rtot1 = %.9f'%R1nm.real.sum()))
+ print(colored_T('Ttot1 = %.9f'%T1nm.real.sum()))
+ print(colored_A('Atot = %.9f'%Q.sum()))
+ print(colored_TOT('TOTAL1 = %.9f'%TOT1))
+ print(colored_i('===> Computed from diffraction efficiencies with normal component'))
+ print(colored_R('Rtot2 = %.9f'%R2nm.real.sum()))
+ print(colored_T('Ttot2 = %.9f'%T2nm.real.sum()))
+ print(colored_A('Atot = %.9f'%Q.sum()))
+ print(colored_TOT('TOTAL2 = %.9f'%TOT2))
+ print(colored_i('===> Computed manually (trapz) from normal component of Poynting vector'))
+ print(colored_i('(expected to be less precise)'))
+ print(colored_R('Rtot3 = %.9f'%R_poy))
+ print(colored_T('Ttot3 = %.9f'%T_poy))
+ print(colored_A('Atot = %.9f'%Q.sum()))
+ print(colored_TOT('TOTAL3 = %.9f'%(Q.sum()+R_poy+T_poy)))
diff --git a/DiffractionGratings/grating3D_runall.sh b/DiffractionGratings/grating3D_runall.sh
index b296018..03c3d6f 100644
--- a/DiffractionGratings/grating3D_runall.sh
+++ b/DiffractionGratings/grating3D_runall.sh
@@ -1,7 +1,8 @@
rm -r res3D*
# for t in bisin checker halfellipsoid hole pyramid torus 2Dlamellar solarcell conv skew
-for t in halfellipsoid torus
+# for t in halfellipsoid torus
+for t in halfellipsoid
do
for f in 0 1
do
@@ -11,7 +12,7 @@ do
done
# for t in bisin checker halfellipsoid hole pyramid torus 2Dlamellar solarcell conv skew
-for t in halfellipsoid torus
+for t in halfellipsoid
do
for f in 0 1
do
@@ -19,4 +20,4 @@ do
echo "____________\ntest case "$t" - FLAG_TOT="$f
python grating3D_postplot.py res3D_$t"_FLAG_TOTAL_"$f
done
-done
\ No newline at end of file
+done
--
GitLab