diff --git a/DiffractionGratings/grating3D.geo b/DiffractionGratings/grating3D.geo
index 8b9384c31d019f943f633daaa05278749dd90844..ea1907ef80bc5e77d7d50cc894917c1311d9e298 100644
--- a/DiffractionGratings/grating3D.geo
+++ b/DiffractionGratings/grating3D.geo
@@ -15,118 +15,131 @@
// You should have received a copy of the GNU General Public License
// along with This program. If not, see <https://www.gnu.org/licenses/>.
-Include "grating3D_data.geo";
+ Include "grating3D_data.geo";
-SetFactory("OpenCASCADE");
-e = 5*nm;
-E = 10*(period_x+period_x/2);
-
-Macro SetPBCs
- BlochXm() = Surface In BoundingBox{.5*(-period_x-dys)-e,-dyc/2-e, PML_bot_hh-e,(-period_x+dys)/2+e, dyc/2+e, PML_top_hh+PML_top+e};
- BlochXp() = Surface In BoundingBox{.5*( period_x-dys)-e,-dyc/2-e, PML_bot_hh-e,( period_x+dys)/2+e, dyc/2+e, PML_top_hh+PML_top+e};
- BlochYm() = Surface In BoundingBox{.5*(-period_x-dys)-e,-dyc/2-e, PML_bot_hh-e,( period_x+dys)/2+e,-dyc/2+e, PML_top_hh+PML_top+e};
- BlochYp() = Surface In BoundingBox{.5*(-period_x-dys)-e, dyc/2-e, PML_bot_hh-e,( period_x+dys)/2+e, dyc/2+e, PML_top_hh+PML_top+e};
- // For k In {1:#BlochXm()-1}
- // Printf("BlochXm surf %g",BlochXm(k));
- // EndFor
- // For k In {1:#BlochYm()-1}
- // Printf("BlochYm surf %g",BlochYm(k));
- // EndFor
- // For k In {1:#BlochXp()-1}
- // Printf("BlochXp surf %g",BlochXp(k));
- // EndFor
- // For k In {1:#BlochYp()-1}
- // Printf("BlochYp surf %g",BlochYp(k));
- // EndFor
- If (tag_geom==6) // bi-sin : BoundingBox does not catch BSpline Surfaces?
- BlochXm()+={17,12};
- BlochXp()+={19,14};
- BlochYm()+={15,20};
- BlochYp()+={13,18};
- EndIf
- Periodic Surface{BlochXp()} = {BlochXm()} Translate{period_x, 0, 0};
- Periodic Surface{BlochYp()} = {BlochYm()} Translate{ dys, dyc, 0};
-Return
-
-If (tag_geom==6)
- // Bi-sinusoidal profile by Spline surface filling
- // not the most general way to implement a freesurface height=f(x,y) and relies on multiple Coherence calls
- // several hacks are used, one of them is to do all this at the very begining
- nsin = 20;
- For i In {0:nsin}
- xx1~{i} = -period_x/2+i*period_x/2/(nsin-1);
- xx2~{i} = i*period_x/2/(nsin-1);
- EndFor
- pp=newp;//pp=pp-1;
- For i In {0:nsin-1}
- yy=-period_y/2;
- Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
- Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
- yy=0;
- Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
- Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
- yy=period_y/2;
- Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
- Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
- yy=-period_y/2;
- Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
- Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
- yy=0;
- Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
- Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
- yy=period_y/2;
- Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
- Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
- EndFor
- For k In {1:11}
- BSpline(newl) = {pp+k:pp+12*nsin:12};
- EndFor
- BSpline(newl) = {pp :pp+12*(nsin-1):12};
- Coherence;
- Box(1) = {-period_x/2,-period_y/2, hh_L_3,period_x,period_y, thick_L_3};
- Delete{Volume{1}; Surface{1,2,3,4} ; Line{13,15,17,19};}
- Curve Loop(7) = {21, 20, -23, -16};
- Curve Loop(8) = {1, 2, -3, -12};
- Surface(7) = {8};
- Curve Loop(10) = {6, 5, -8, -3};
- Surface(8) = {10};
- Curve Loop(12) = {9, 10, -11, -8};
- Surface(9) = {12};
- Curve Loop(14) = {2, 9, -4, -7};
- Surface(10) = {14};
- Line(25) = {239, 229};
- Line(26) = {229, 238};
- Line(27) = {243, 235};
- Line(28) = {235, 242};
- Line(29) = {237, 244};
- Line(30) = {233, 240};
- Line(31) = {237, 245};
- Line(32) = {241, 233};
- Curve Loop(16) = {25, 12, 6, -27, -21};
- Surface(11) = {16};
- Curve Loop(18) = {20, -31, -11, -5, -27};
- Surface(12) = {18};
- Curve Loop(20) = {29, -18, -28, 5, 11};
- Surface(13) = {20};
- Curve Loop(22) = {10, 29, -24, -30, 4};
- Surface(14) = {22};
- Curve Loop(24) = {30, -14, -26, 1, 7};
- Surface(15) = {24};
- Curve Loop(26) = {25, 1, 7, -32, -16};
- Surface(16) = {26};
- Curve Loop(28) = {26, 22, -28, -6, -12};
- Surface(17) = {28};
- Curve Loop(30) = {31, -23, 32, 4, 10};
- Surface(18) = {30};
- Surface Loop(3) = {5, 16, 18, 12, 11, 9, 8, 7, 10};
- Volume(1) = {3};
- Surface Loop(2) = {6, 15, 14, 10, 7, 8, 9, 13, 17};
- Volume(2) = {2};
-EndIf
-
-For k In {7:0:-1}
+ SetFactory("OpenCASCADE");
+ e = 5*nm;
+ E = 10*(period_x+period_x/2);
+
+ Macro SetPBCs
+ BlochXm() = Surface In BoundingBox{.5*(-period_x-dys)-e,-dyc/2-e, PML_bot_hh-e,(-period_x+dys)/2+e, dyc/2+e, PML_top_hh+PML_top+e};
+ BlochXp() = Surface In BoundingBox{.5*( period_x-dys)-e,-dyc/2-e, PML_bot_hh-e,( period_x+dys)/2+e, dyc/2+e, PML_top_hh+PML_top+e};
+ BlochYm() = Surface In BoundingBox{.5*(-period_x-dys)-e,-dyc/2-e, PML_bot_hh-e,( period_x+dys)/2+e,-dyc/2+e, PML_top_hh+PML_top+e};
+ BlochYp() = Surface In BoundingBox{.5*(-period_x-dys)-e, dyc/2-e, PML_bot_hh-e,( period_x+dys)/2+e, dyc/2+e, PML_top_hh+PML_top+e};
+ // For k In {1:#BlochXm()-1}
+ // Printf("BlochXm surf %g",BlochXm(k));
+ // EndFor
+ // For k In {1:#BlochYm()-1}
+ // Printf("BlochYm surf %g",BlochYm(k));
+ // EndFor
+ // For k In {1:#BlochXp()-1}
+ // Printf("BlochXp surf %g",BlochXp(k));
+ // EndFor
+ // For k In {1:#BlochYp()-1}
+ // Printf("BlochYp surf %g",BlochYp(k));
+ // EndFor
+ If (tag_geom==6) // bi-sin : BoundingBox does not catch BSpline Surfaces?
+ BlochXm()+={17,12};
+ BlochXp()+={19,14};
+ BlochYm()+={15,20};
+ BlochYp()+={13,18};
+ EndIf
+ Periodic Surface{BlochXp()} = {BlochXm()} Translate{period_x, 0, 0};
+ Periodic Surface{BlochYp()} = {BlochYm()} Translate{ dys, dyc, 0};
+ Return
+
If (tag_geom==6)
- If (k!=3)
+ // Bi-sinusoidal profile by Spline surface filling
+ // not the most general way to implement a freesurface height=f(x,y) and relies on multiple Coherence calls
+ // several hacks are used, one of them is to do all this at the very begining
+ nsin = 20;
+ For i In {0:nsin}
+ xx1~{i} = -period_x/2+i*period_x/2/(nsin-1);
+ xx2~{i} = i*period_x/2/(nsin-1);
+ EndFor
+ pp=newp;//pp=pp-1;
+ For i In {0:nsin-1}
+ yy=-period_y/2;
+ Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+ Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
+ yy=0;
+ Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+ Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
+ yy=period_y/2;
+ Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+ Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
+ yy=-period_y/2;
+ Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+ Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
+ yy=0;
+ Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+ Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
+ yy=period_y/2;
+ Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+ Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
+ EndFor
+ For k In {1:11}
+ BSpline(newl) = {pp+k:pp+12*nsin:12};
+ EndFor
+ BSpline(newl) = {pp :pp+12*(nsin-1):12};
+ Coherence;
+ Box(1) = {-period_x/2,-period_y/2, hh_L_3,period_x,period_y, thick_L_3};
+ Delete{Volume{1}; Surface{1,2,3,4} ; Line{13,15,17,19};}
+ Curve Loop(7) = {21, 20, -23, -16};
+ Curve Loop(8) = {1, 2, -3, -12};
+ Surface(7) = {8};
+ Curve Loop(10) = {6, 5, -8, -3};
+ Surface(8) = {10};
+ Curve Loop(12) = {9, 10, -11, -8};
+ Surface(9) = {12};
+ Curve Loop(14) = {2, 9, -4, -7};
+ Surface(10) = {14};
+ Line(25) = {239, 229};
+ Line(26) = {229, 238};
+ Line(27) = {243, 235};
+ Line(28) = {235, 242};
+ Line(29) = {237, 244};
+ Line(30) = {233, 240};
+ Line(31) = {237, 245};
+ Line(32) = {241, 233};
+ Curve Loop(16) = {25, 12, 6, -27, -21};
+ Surface(11) = {16};
+ Curve Loop(18) = {20, -31, -11, -5, -27};
+ Surface(12) = {18};
+ Curve Loop(20) = {29, -18, -28, 5, 11};
+ Surface(13) = {20};
+ Curve Loop(22) = {10, 29, -24, -30, 4};
+ Surface(14) = {22};
+ Curve Loop(24) = {30, -14, -26, 1, 7};
+ Surface(15) = {24};
+ Curve Loop(26) = {25, 1, 7, -32, -16};
+ Surface(16) = {26};
+ Curve Loop(28) = {26, 22, -28, -6, -12};
+ Surface(17) = {28};
+ Curve Loop(30) = {31, -23, 32, 4, 10};
+ Surface(18) = {30};
+ Surface Loop(3) = {5, 16, 18, 12, 11, 9, 8, 7, 10};
+ Volume(1) = {3};
+ Surface Loop(2) = {6, 15, 14, 10, 7, 8, 9, 13, 17};
+ Volume(2) = {2};
+ EndIf
+
+ For k In {7:0:-1}
+ If (tag_geom==6)
+ If (k!=3)
+ p=newp; l=newl; ll=newll; s=news;
+ Point(p) = {0.5*(-period_x-dys), -dyc/2, hh_L~{k}};
+ Point(p+1) = {0.5*( period_x-dys), -dyc/2, hh_L~{k}};
+ Point(p+2) = {0.5*( period_x+dys), dyc/2, hh_L~{k}};
+ Point(p+3) = {0.5*(-period_x+dys), dyc/2, hh_L~{k}};
+ Line(l)={p,p+1};Line(l+1)={p+1,p+2};Line(l+2)={p+2,p+3};Line(l+3)={p+3,p};
+ Curve Loop(ll) = {l,l+1,l+2,l+3};
+ Surface(s) = {ll};
+ Extrude {0, 0, thick_L~{k}} {Surface{s};}
+ Else
+ v=newv;
+ EndIf
+ Else
p=newp; l=newl; ll=newll; s=news;
Point(p) = {0.5*(-period_x-dys), -dyc/2, hh_L~{k}};
Point(p+1) = {0.5*( period_x-dys), -dyc/2, hh_L~{k}};
@@ -136,173 +149,161 @@ For k In {7:0:-1}
Curve Loop(ll) = {l,l+1,l+2,l+3};
Surface(s) = {ll};
Extrude {0, 0, thick_L~{k}} {Surface{s};}
- Else
- v=newv;
EndIf
+ EndFor
+ // Box(1) = {-period_x/2,-period_y/2, PML_bot_hh,period_x,period_y, PML_bot };
+ // Box(2) = {-period_x/2,-period_y/2, hh_L_6,period_x,period_y, thick_L_6};
+ // Box(3) = {-period_x/2,-period_y/2, hh_L_5,period_x,period_y, thick_L_5};
+ // Box(4) = {-period_x/2,-period_y/2, hh_L_4,period_x,period_y, thick_L_4};
+ // If (tag_geom!=6)
+ // Box(5) = {-period_x/2,-period_y/2, hh_L_3,period_x,period_y, thick_L_3};
+ // EndIf
+ // Box(6) = {-period_x/2,-period_y/2, hh_L_2,period_x,period_y, thick_L_2};
+ // Box(7) = {-period_x/2,-period_y/2, hh_L_1,period_x,period_y, thick_L_1};
+ // Box(8) = {-period_x/2,-period_y/2, PML_top_hh,period_x,period_y, PML_top };
+
+ If (tag_geom==1)
+ p=newp;
+ Point(p) = {0,0,hh_L_3+rz};
+ Line(97) = {29, 65};
+ Line(98) = {65, 32};
+ Line(99) = {65, 30};
+ Line(100) = {65, 31};
+ Curve Loop(92) = {97, 99, -43};
+ Plane Surface(91) = {92};
+ Curve Loop(93) = {99, 45, -100};
+ Plane Surface(92) = {93};
+ Curve Loop(94) = {100, 47, -98};
+ Plane Surface(93) = {94};
+ Curve Loop(95) = {97, 98, 48};
+ Plane Surface(94) = {95};
+ Surface Loop(9) = {91, 94, 93, 92, 42};
+ Volume(9) = {9};
+ EndIf
+
+ If (tag_geom==2)
+ Cylinder(9) = {0,0,hh_L_3,0,0,thick_L_3,rx};
+ EndIf
+
+ If (tag_geom==3)
+ Lx = 210*nm;
+ Ly = 232*nm;
+ Box(9) = {-Lx/2 , -Ly/2 , hh_L_3 , Lx, ry, rz};
+ Box(10) = {-Lx/2 , Ly/2-ry , hh_L_3 , Lx, ry, rz};
+ Box(11) = {-Lx/2 , -Ly/2 , hh_L_3 , rx, Ly, rz};
+ scat() = BooleanUnion{ Volume{9}; Delete; }{ Volume{10,11}; Delete;} ;
+ r_fillet = 10*nm;
+ fscat() = Abs(Boundary{ Volume{scat()} ; });
+ escat() = Unique(Abs(Boundary{ Surface{fscat()}; }));
+ // Fillet{scat()}{escat()}{r_fillet}
+ EndIf
+
+ If (tag_geom==4)
+ Sphere(9) = {0,0,hh_L_3,rx};
+ Dilate { { 0,0,hh_L_3 }, { 1, ry/rx, rz/rx } } { Volume{9}; }
+ BooleanDifference{ Volume{9}; Delete; }{ Volume{4}; }
+ EndIf
+
+ If (tag_geom==5)
+ Box(9) = {-rx/2,-ry/2, hh_L_2-rz, rx, ry, rz};
+ Rotate {{0, 0, 1}, {0,0,0}, Pi/4} {Volume{9};}
+ EndIf
+
+ If (tag_geom==7)
+ Box(9) = {-period_x/2,-ry/2, hh_L_3, period_x, ry, rz};
+ EndIf
+
+
+ Coherence;
+ Call SetPBCs;
+ If (tag_geom==6)
+ Physical Volume("PMLBOT" ,1) = {3};
+ Physical Volume("LAYER_L6_SUBS" ,2) = {4};
+ Physical Volume("LAYER_L5" ,3) = {5};
+ Physical Volume("LAYER_L4" ,4) = {6};
+ Physical Volume("LAYER_L3" ,5) = {2};
+ Physical Volume("LAYER_L2" ,6) = {7};
+ Physical Volume("LAYER_L1_SUPER",7) = {8};
+ Physical Volume("PMLTOP" ,8) = {9};
+ Physical Volume("SCAT" ,9) = {1};
Else
- p=newp; l=newl; ll=newll; s=news;
- Point(p) = {0.5*(-period_x-dys), -dyc/2, hh_L~{k}};
- Point(p+1) = {0.5*( period_x-dys), -dyc/2, hh_L~{k}};
- Point(p+2) = {0.5*( period_x+dys), dyc/2, hh_L~{k}};
- Point(p+3) = {0.5*(-period_x+dys), dyc/2, hh_L~{k}};
- Line(l)={p,p+1};Line(l+1)={p+1,p+2};Line(l+2)={p+2,p+3};Line(l+3)={p+3,p};
- Curve Loop(ll) = {l,l+1,l+2,l+3};
- Surface(s) = {ll};
- Extrude {0, 0, thick_L~{k}} {Surface{s};}
+ Physical Volume("PMLBOT" ,1) = {1};
+ Physical Volume("LAYER_L6_SUBS" ,2) = {2};
+ Physical Volume("LAYER_L5" ,3) = {3};
+ Physical Volume("LAYER_L4" ,4) = {4};
+ Physical Volume("LAYER_L3" ,5) = {10};
+ Physical Volume("LAYER_L2" ,6) = {6};
+ Physical Volume("LAYER_L1_SUPER",7) = {7};
+ Physical Volume("PMLTOP" ,8) = {8};
+ Physical Volume("SCAT" ,9) = {9};
EndIf
-EndFor
-// Box(1) = {-period_x/2,-period_y/2, PML_bot_hh,period_x,period_y, PML_bot };
-// Box(2) = {-period_x/2,-period_y/2, hh_L_6,period_x,period_y, thick_L_6};
-// Box(3) = {-period_x/2,-period_y/2, hh_L_5,period_x,period_y, thick_L_5};
-// Box(4) = {-period_x/2,-period_y/2, hh_L_4,period_x,period_y, thick_L_4};
-// If (tag_geom!=6)
-// Box(5) = {-period_x/2,-period_y/2, hh_L_3,period_x,period_y, thick_L_3};
-// EndIf
-// Box(6) = {-period_x/2,-period_y/2, hh_L_2,period_x,period_y, thick_L_2};
-// Box(7) = {-period_x/2,-period_y/2, hh_L_1,period_x,period_y, thick_L_1};
-// Box(8) = {-period_x/2,-period_y/2, PML_top_hh,period_x,period_y, PML_top };
-
-If (tag_geom==1)
- p=newp;
- Point(p) = {0,0,hh_L_3+rz};
- Line(97) = {29, 65};
- Line(98) = {65, 32};
- Line(99) = {65, 30};
- Line(100) = {65, 31};
- Curve Loop(92) = {97, 99, -43};
- Plane Surface(91) = {92};
- Curve Loop(93) = {99, 45, -100};
- Plane Surface(92) = {93};
- Curve Loop(94) = {100, 47, -98};
- Plane Surface(93) = {94};
- Curve Loop(95) = {97, 98, 48};
- Plane Surface(94) = {95};
- Surface Loop(9) = {91, 94, 93, 92, 42};
- Volume(9) = {9};
-EndIf
-
-If (tag_geom==2)
- Cylinder(9) = {0,0,hh_L_3,0,0,thick_L_3,rx};
-EndIf
-
-If (tag_geom==3)
- Lx = 210*nm;
- Ly = 232*nm;
- Box(9) = {-Lx/2 , -Ly/2 , hh_L_3 , Lx, ry, rz};
- Box(10) = {-Lx/2 , Ly/2-ry , hh_L_3 , Lx, ry, rz};
- Box(11) = {-Lx/2 , -Ly/2 , hh_L_3 , rx, Ly, rz};
- scat() = BooleanUnion{ Volume{9}; Delete; }{ Volume{10,11}; Delete;} ;
- r_fillet = 10*nm;
- fscat() = Abs(Boundary{ Volume{scat()} ; });
- escat() = Unique(Abs(Boundary{ Surface{fscat()}; }));
- // Fillet{scat()}{escat()}{r_fillet}
-EndIf
-
-If (tag_geom==4)
- Sphere(9) = {0,0,hh_L_3,rx};
- Dilate { { 0,0,hh_L_3 }, { 1, ry/rx, rz/rx } } { Volume{9}; }
- BooleanDifference{ Volume{9}; Delete; }{ Volume{4}; }
-EndIf
-
-If (tag_geom==5)
- Box(9) = {-rx/2,-ry/2, hh_L_2-rz, rx, ry, rz};
- Rotate {{0, 0, 1}, {0,0,0}, Pi/4} {Volume{9};}
-EndIf
-
-If (tag_geom==7)
- Box(9) = {-period_x/2,-ry/2, hh_L_3, period_x, ry, rz};
-EndIf
-
-
-Coherence;
-Call SetPBCs;
-If (tag_geom==6)
- Physical Volume("PMLBOT" ,1) = {3};
- Physical Volume("LAYER_L6_SUBS" ,2) = {4};
- Physical Volume("LAYER_L5" ,3) = {5};
- Physical Volume("LAYER_L4" ,4) = {6};
- Physical Volume("LAYER_L3" ,5) = {2};
- Physical Volume("LAYER_L2" ,6) = {7};
- Physical Volume("LAYER_L1_SUPER",7) = {8};
- Physical Volume("PMLTOP" ,8) = {9};
- Physical Volume("SCAT" ,9) = {1};
-Else
- Physical Volume("PMLBOT" ,1) = {1};
- Physical Volume("LAYER_L6_SUBS" ,2) = {2};
- Physical Volume("LAYER_L5" ,3) = {3};
- Physical Volume("LAYER_L4" ,4) = {4};
- Physical Volume("LAYER_L3" ,5) = {10};
- Physical Volume("LAYER_L2" ,6) = {6};
- Physical Volume("LAYER_L1_SUPER",7) = {7};
- Physical Volume("PMLTOP" ,8) = {8};
- Physical Volume("SCAT" ,9) = {9};
-EndIf
-Physical Surface("BXM" ,101 ) = BlochXm();
-Physical Surface("BXP" ,102 ) = BlochXp();
-Physical Surface("BYM" ,201 ) = BlochYm();
-Physical Surface("BYP" ,202 ) = BlochYp();
-Physical Surface("STOP",301 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, PML_top_hh-e,period_x/2+E, period_y/2+E, PML_top_hh+e};
-Physical Surface("SBOT",302 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, hh_L_6-e,period_x/2+E , period_y/2+E, hh_L_6+e};
-Physical Surface("SPMLTOP",401 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, PML_top_hh+PML_top-e,period_x/2+E, period_y/2+E, PML_top_hh+PML_top+e};
-Physical Surface("SPMLBOT",402 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, PML_bot_hh-e,period_x/2+E, period_y/2+E, PML_bot_hh+e};
-Physical Surface("SVIS",500 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, hh_L_3-e,period_x/2+E , period_y/2+E, hh_L_3+e};
-
-pts_PMLBOT() = PointsOf{ Physical Volume{1}; };
-pts_LAYER_L6() = PointsOf{ Physical Volume{2}; };
-pts_LAYER_L5() = PointsOf{ Physical Volume{3}; };
-pts_LAYER_L4() = PointsOf{ Physical Volume{4}; };
-pts_LAYER_L3() = PointsOf{ Physical Volume{5}; };
-pts_LAYER_L2() = PointsOf{ Physical Volume{6}; };
-pts_LAYER_L1() = PointsOf{ Physical Volume{7}; };
-pts_PMLTOP() = PointsOf{ Physical Volume{8}; };
-pts_ROD() = PointsOf{ Physical Volume{9}; };
-
-// if test_case==conv, we want to remesh the scatterer as well
-// the following enforces to update lc_scat when paramaille changes
-If (tag_geom==4)
- lc_scat = lambda_m/(paramaille*5);
-EndIf
-
-lc_PML = lambda_m/(paramaille*.5);
-list_lc(0) = lc_PML;
-For k In {1:6}
- n_L~{k} = Sqrt(Abs(eps_re_L~{k}));
- lc_L~{k} = lambda_m/(paramaille*n_L~{k})/refine_mesh_L~{k};
- list_lc(7-k) = lc_L~{k};
-EndFor
-list_lc(7) = lc_PML;
-list_lc(8) = lc_scat;
-
-// This helps meshing: The default behavior of the PointsOf technique
-// overides points belonging to several domains (by order of call of Characteristic Length)
-If (tag_geom==7)
- meshing_sequence() = {1,8,2,3,5,6,7,4,9};
-Else
- meshing_sequence() = {1,8,2,3,4,6,7,5,9};
-EndIf
-
-// Start with coarsest
-Characteristic Length{:} = lc_PML;
-
-For k In {0:8}
- which_dom = meshing_sequence(k);
- Characteristic Length{PointsOf{Physical Volume{which_dom};}} = list_lc(which_dom-1);
-EndFor
-
-
-If (tag_geom==3) // Split U weird otherwise
- Mesh.Algorithm = 6;
-EndIf
-
-
-// Mesh.SurfaceEdges = 0;
-Mesh.VolumeEdges = 0;
-Mesh.ElementOrder = og;
-If (og==2)
- Mesh.HighOrderOptimize = 1;
-EndIf
-If (og==1)
- Mesh.Optimize = 1;
-EndIf
-
-// Solver.SocketName=Sprintf("socktest%g",socketid);
+ Physical Surface("BXM" ,101 ) = BlochXm();
+ Physical Surface("BXP" ,102 ) = BlochXp();
+ Physical Surface("BYM" ,201 ) = BlochYm();
+ Physical Surface("BYP" ,202 ) = BlochYp();
+ Physical Surface("STOP",301 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, PML_top_hh-e,period_x/2+E, period_y/2+E, PML_top_hh+e};
+ Physical Surface("SBOT",302 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, hh_L_6-e,period_x/2+E , period_y/2+E, hh_L_6+e};
+ Physical Surface("SPMLTOP",401 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, PML_top_hh+PML_top-e,period_x/2+E, period_y/2+E, PML_top_hh+PML_top+e};
+ Physical Surface("SPMLBOT",402 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, PML_bot_hh-e,period_x/2+E, period_y/2+E, PML_bot_hh+e};
+ Physical Surface("SVIS",500 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, hh_L_3-e,period_x/2+E , period_y/2+E, hh_L_3+e};
+
+ pts_PMLBOT() = PointsOf{ Physical Volume{1}; };
+ pts_LAYER_L6() = PointsOf{ Physical Volume{2}; };
+ pts_LAYER_L5() = PointsOf{ Physical Volume{3}; };
+ pts_LAYER_L4() = PointsOf{ Physical Volume{4}; };
+ pts_LAYER_L3() = PointsOf{ Physical Volume{5}; };
+ pts_LAYER_L2() = PointsOf{ Physical Volume{6}; };
+ pts_LAYER_L1() = PointsOf{ Physical Volume{7}; };
+ pts_PMLTOP() = PointsOf{ Physical Volume{8}; };
+ pts_ROD() = PointsOf{ Physical Volume{9}; };
+
+ // if test_case==conv, we want to remesh the scatterer as well
+ // the following enforces to update lc_scat when paramaille changes
+ If (tag_geom==4)
+ lc_scat = lambda_m/(paramaille*5);
+ EndIf
+
+ lc_PML = lambda_m/(paramaille*.5);
+ list_lc(0) = lc_PML;
+ For k In {1:6}
+ n_L~{k} = Sqrt(Abs(eps_re_L~{k}));
+ lc_L~{k} = lambda_m/(paramaille*n_L~{k})/refine_mesh_L~{k};
+ list_lc(7-k) = lc_L~{k};
+ EndFor
+ list_lc(7) = lc_PML;
+ list_lc(8) = lc_scat;
+
+ // This helps meshing: The default behavior of the PointsOf technique
+ // overides points belonging to several domains (by order of call of Characteristic Length)
+ If (tag_geom==7)
+ meshing_sequence() = {1,8,2,3,5,6,7,4,9};
+ Else
+ meshing_sequence() = {1,8,2,3,4,6,7,5,9};
+ EndIf
+
+ // Start with coarsest
+ Characteristic Length{:} = lc_PML;
+
+ For k In {0:8}
+ which_dom = meshing_sequence(k);
+ Characteristic Length{PointsOf{Physical Volume{which_dom};}} = list_lc(which_dom-1);
+ EndFor
+
+
+ // If (tag_geom==3) // Split U weird otherwise
+ Mesh.Algorithm = 1;
+ // EndIf
+
+
+ // Mesh.SurfaceEdges = 0;
+ Mesh.VolumeEdges = 0;
+ Mesh.ElementOrder = og;
+ If (og==2)
+ Mesh.HighOrderOptimize = 1;
+ EndIf
+ If (og==1)
+ Mesh.Optimize = 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 0348ee7954691006089154c2aae08569146b65d8..14bf789b86c7df77d6f2ddedd1dca9eda4119433 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -15,541 +15,561 @@
// You should have received a copy of the GNU General Public License
// along with This program. If not, see <https://www.gnu.org/licenses/>.
-Include "grating3D_data.geo";
-Include "grating3D_materials.pro";
-
-myDir = "res3D/";
-
-Group {
- // SubDomains
- PMLbot = Region[1];
- L_6_temp = Region[2];
- L_5 = Region[3];
- L_4 = Region[4];
- L_3 = Region[5];
- L_2 = Region[6];
- L_1_temp = Region[7];
- PMLtop = Region[8];
- Scat = Region[9];
-
- // Boundaries
- SurfBlochXm = Region[101];
- SurfBlochXp = Region[102];
- SurfBlochYm = Region[201];
- SurfBlochYp = Region[202];
- SurfIntTop = Region[301];
- SurfIntBot = Region[302];
-
-
- SurfDirichlet = Region[{401,402}];
- SurfBloch = Region[{SurfBlochXm,SurfBlochXp,SurfBlochYm,SurfBlochYp}];
-
- If (FlagLinkFacets==1)
- SurfExcludeFacets = Region[{}];
- Else
- SurfExcludeFacets = Region[{SurfBloch}];
- EndIf
-
- 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}];
-
- // get normal comp of E field on integ surfaces
- Gama = Region[{SurfIntBot,SurfIntTop}];
- Tr = ElementsOf[Omega_plot, ConnectedTo Gama];
-}
-
-
-
-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;
- cel = 1/Sqrt[ep0 * mu0];
- om0 = 2*Pi*cel/lambda0;
- k0 = 2.*Pi/lambda0;
- Ae = 1;
- Pinc = 0.5*Ae^2*Sqrt[eps_re_L_1*ep0/mu0] * Cos[theta0];
-
- // permittivities
- For i In {1:6}
- If (flag_mat~{i}==0)
- epsr[L~{i}] = Complex[eps_re_L~{i} , eps_im_L~{i}] * TensorDiag[1,1,1];
- If (i==1)
- epsr1[] = Complex[eps_re_L~{i} , eps_im_L~{i}];
- EndIf
- If (i==6)
- epsr2[] = Complex[eps_re_L~{i} , eps_im_L~{i}];
- EndIf
+ Include "grating3D_data.geo";
+ Include "grating3D_materials.pro";
+
+ myDir = "res3D/";
+
+ Group {
+ // SubDomains
+ PMLbot = Region[1];
+ L_6 = Region[2];
+ L_5 = Region[3];
+ L_4 = Region[4];
+ L_3 = Region[5];
+ L_2 = Region[6];
+ L_1 = Region[7];
+ PMLtop = Region[8];
+ Scat = Region[9];
+
+ // Boundaries
+ SurfBlochXm = Region[101];
+ SurfBlochXp = Region[102];
+ SurfBlochYm = Region[201];
+ SurfBlochYp = Region[202];
+ SurfIntTop = Region[301];
+ SurfIntBot = Region[302];
+
+
+ SurfDirichlet = Region[{401,402}];
+ SurfBloch = Region[{SurfBlochXm,SurfBlochXp,SurfBlochYm,SurfBlochYp}];
+
+ If (FlagLinkFacets==1)
+ SurfExcludeFacets = Region[{}];
Else
- For j In {1:nb_available_materials}
- If(flag_mat~{i}==j)
- epsr[L~{i}] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]] * TensorDiag[1,1,1];
- If (i==1)
- epsr1[] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]];
- EndIf
- If (i==6)
- epsr2[] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]];
- EndIf
+ SurfExcludeFacets = Region[{SurfBloch}];
+ EndIf
+
+ // 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}];
+
+ // get normal comp of E field on integ surfaces
+ Gama = Region[{SurfIntBot,SurfIntTop}];
+ Tr = ElementsOf[Omega_plot, ConnectedTo Gama];
+ }
+
+
+
+ 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;
+ cel = 1/Sqrt[ep0 * mu0];
+ om0 = 2*Pi*cel/lambda0;
+ k0 = 2.*Pi/lambda0;
+ Ae = 1;
+ Pinc = 0.5*Ae^2*Sqrt[eps_re_L_1*ep0/mu0] * Cos[theta0];
+
+ // permittivities
+ For i In {1:6}
+ If (flag_mat~{i}==0)
+ epsr[L~{i}] = Complex[eps_re_L~{i} , eps_im_L~{i}] * TensorDiag[1,1,1];
+ If (i==1)
+ epsr1[] = Complex[eps_re_L~{i} , eps_im_L~{i}];
+ EndIf
+ If (i==6)
+ epsr2[] = Complex[eps_re_L~{i} , eps_im_L~{i}];
EndIf
+ Else
+ For j In {1:nb_available_materials}
+ If(flag_mat~{i}==j)
+ epsr[L~{i}] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]] * TensorDiag[1,1,1];
+ If (i==1)
+ epsr1[] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]];
+ EndIf
+ If (i==6)
+ epsr2[] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]];
+ EndIf
+ EndIf
+ EndFor
+ EndIf
+ EndFor
+
+ If (flag_mat_scat==0)
+ epsr[Scat] = Complex[eps_re_Scat , eps_im_Scat] * TensorDiag[1,1,1];
+ Else
+ For j In {flag_mat_scat:flag_mat_scat}
+ epsr[Scat] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]] * TensorDiag[1,1,1];
EndFor
EndIf
- EndFor
-
- If (flag_mat_scat==0)
- epsr[Scat] = Complex[eps_re_Scat , eps_im_Scat] * TensorDiag[1,1,1];
- Else
- For j In {flag_mat_scat:flag_mat_scat}
- epsr[Scat] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]] * TensorDiag[1,1,1];
+
+ For i In {1:6}
+ mur[L~{i}] = TensorDiag[1,1,1];
EndFor
- EndIf
-
- For i In {1:6}
- mur[L~{i}] = TensorDiag[1,1,1];
- EndFor
- mur[Scat] = TensorDiag[1,1,1];
-
- ////// PMLS
- a_pml = 1.;
- b_pml = 1.;
- // bermu
- n1[] = Sqrt[epsr1[]];
- n2[] = Sqrt[epsr2[]];
- k1norm[] = k0*n1[];
- k2norm[] = k0*n2[];
-
- Zmax = PML_top_hh;
- Zmin = hh_L_6;
- Damp_pml_top[] = 1/(Zmax + PML_top - Fabs[Z[]]) - 1/(PML_top);
- Damp_pml_bot[] = 1/(Zmin + PML_top - Fabs[Z[]]) - 1/(PML_bot);
- Sigma_top[] = 0.5*(Damp_pml_top[] + Fabs[Damp_pml_top[]]);
- Sigma_bot[] = 0.5*(Damp_pml_bot[] + Fabs[Damp_pml_bot[]]);
-
- If (PML_TYPE==0)
- sz[] = Complex[a_pml,b_pml];
- Else
- sz[PMLtop] = Complex[1,Damp_pml_top[]/k1norm[]];
- sz[PMLbot] = Complex[1,Damp_pml_bot[]/k2norm[]];
- EndIf
- sx = 1.;
- sy = 1.;
-
- epsr[PMLtop] = Re[epsr1[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
- mur[PMLtop] = TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
- epsr[PMLbot] = Re[epsr2[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
- mur[PMLbot] = TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
-
- // epsr[PMLtop] = Re[epsr1[]]*TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
- // mur[PMLtop] = TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
- // epsr[PMLbot] = Re[epsr2[]]*TensorDiag[sz_bermubot[]*sy/sx,sx*sz_bermubot[]/sy,sx*sy/sz_bermubot[]];
- // mur[PMLbot] = TensorDiag[sz_bermubot[]*sy/sx,sx*sz_bermubot[]/sy,sx*sy/sz_bermubot[]];
-
- // epsr[PMLbot] = epsr2[];
- // mur[PMLbot] = TensorDiag[1,1,1];
-
- epsr_annex[PMLbot] = epsr[];
- epsr_annex[PMLtop] = epsr[];
- epsr_annex[Omega_source] = epsr1[] * TensorDiag[1,1,1];
- epsr_annex[L_1] = epsr[];
- epsr_annex[L_6] = epsr[];
-
- //// 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];
- ppol_r[] = (k1r[]/Norm[k1r[]]) /\ spol[];
- ppol_t[] = (k2[] /Norm[k2[]] ) /\ spol[];
-
- AmplEis[] = spol[];
- AmplErs[] = rs[]*spol[];
- AmplEts[] = ts[]*spol[];
- AmplHis[] = Sqrt[ep0*epsr1[]/mu0] *spol[];
- AmplHrs[] = Sqrt[ep0*epsr1[]/mu0]*rp[]*spol[];
- AmplHts[] = Sqrt[ep0*epsr1[]/mu0]*tp[]*spol[];
-
- Eis[] = AmplEis[] * Exp[I[]*k1[] *XYZ[]];
- Ers[] = AmplErs[] * Exp[I[]*k1r[]*XYZ[]];
- Ets[] = AmplEts[] * Exp[I[]*k2[] *XYZ[]];
- His[] = AmplHis[] * Exp[I[]*k1[] *XYZ[]];
- Hrs[] = AmplHrs[] * Exp[I[]*k1r[]*XYZ[]];
- Hts[] = AmplHts[] * Exp[I[]*k2[] *XYZ[]];
- Eip[] = -1/(om0*ep0*epsr1[]) * k1[] /\ His[];
- Erp[] = -1/(om0*ep0*epsr1[]) * k1r[] /\ Hrs[];
- Etp[] = -1/(om0*ep0*epsr2[]) * k2[] /\ Hts[];
-
- Ei[] = Ae*(Cos[psi0]*Eip[]-Sin[psi0]*Eis[]);
- Er[] = Ae*(Cos[psi0]*Erp[]-Sin[psi0]*Ers[]);
- Et[] = Ae*(Cos[psi0]*Etp[]-Sin[psi0]*Ets[]);
- Hi[] = 1/(om0*mu0*mur[]) * k1[] /\ Ei[];
- Hr[] = 1/(om0*mu0*mur[]) * k1r[] /\ Er[];
- Ht[] = 1/(om0*mu0*mur[]) * k2[] /\ Et[];
-
- E1[Omega_super] = Ei[]+Er[];
- E1[Omega_subs] = Et[];
- E1d[Omega_super] = Er[];
- E1d[Omega_subs] = Et[];
- H1[Omega_super] = Hi[]+Hr[];
- H1[Omega_subs] = Ht[];
- H1d[Omega_super] = Hr[];
- H1d[Omega_subs] = Ht[];
-
- source_vol_scat[] = (om0/cel)^2*(epsr[]-epsr_annex[])*E1[];
- source_surf_tot[] = -2*I[]*om0*mu0*Vector[0,0,1] /\ Hi[];
- // 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];
- dephY[] = Exp[I[]*sky1[]*period_y];
-
- // Fourier coefficients variables
- Nb_ordre = 2*Nmax+1;
- For i In {0:Nb_ordre-1}
- For j In {0:Nb_ordre-1}
- 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[])];
+ mur[Scat] = TensorDiag[1,1,1];
+
+ ////// PMLS
+ a_pml = 1.;
+ b_pml = 1.;
+ // bermu
+ n1[] = Sqrt[epsr1[]];
+ n2[] = Sqrt[epsr2[]];
+ k1norm[] = k0*n1[];
+ k2norm[] = k0*n2[];
+
+ Zmax = PML_top_hh;
+ Zmin = hh_L_6;
+ Damp_pml_top[] = 1/(Zmax + PML_top - Fabs[Z[]]) - 1/(PML_top);
+ Damp_pml_bot[] = 1/(Zmin + PML_top - Fabs[Z[]]) - 1/(PML_bot);
+ Sigma_top[] = 0.5*(Damp_pml_top[] + Fabs[Damp_pml_top[]]);
+ Sigma_bot[] = 0.5*(Damp_pml_bot[] + Fabs[Damp_pml_bot[]]);
+
+ If (PML_TYPE==0)
+ sz[] = Complex[a_pml,b_pml];
+ Else
+ sz[PMLtop] = Complex[1,Damp_pml_top[]/k1norm[]];
+ sz[PMLbot] = Complex[1,Damp_pml_bot[]/k2norm[]];
+ EndIf
+ sx = 1.;
+ sy = 1.;
+
+ epsr[PMLtop] = Re[epsr1[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
+ mur[PMLtop] = TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
+ epsr[PMLbot] = Re[epsr2[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
+ mur[PMLbot] = TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
+
+ // epsr[PMLtop] = Re[epsr1[]]*TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
+ // mur[PMLtop] = TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
+ // epsr[PMLbot] = Re[epsr2[]]*TensorDiag[sz_bermubot[]*sy/sx,sx*sz_bermubot[]/sy,sx*sy/sz_bermubot[]];
+ // mur[PMLbot] = TensorDiag[sz_bermubot[]*sy/sx,sx*sz_bermubot[]/sy,sx*sy/sz_bermubot[]];
+
+ // epsr[PMLbot] = epsr2[];
+ // mur[PMLbot] = TensorDiag[1,1,1];
+
+ epsr_annex[PMLbot] = epsr[];
+ epsr_annex[PMLtop] = epsr[];
+ epsr_annex[Omega_source] = epsr1[] * TensorDiag[1,1,1];
+ epsr_annex[L_1] = epsr[];
+ epsr_annex[L_6] = epsr[];
+
+ //// 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];
+ ppol_r[] = (k1r[]/Norm[k1r[]]) /\ spol[];
+ ppol_t[] = (k2[] /Norm[k2[]] ) /\ spol[];
+
+ AmplEis[] = spol[];
+ AmplErs[] = rs[]*spol[];
+ AmplEts[] = ts[]*spol[];
+ AmplHis[] = Sqrt[ep0*epsr1[]/mu0] *spol[];
+ AmplHrs[] = Sqrt[ep0*epsr1[]/mu0]*rp[]*spol[];
+ AmplHts[] = Sqrt[ep0*epsr1[]/mu0]*tp[]*spol[];
+
+ Eis[] = AmplEis[] * Exp[I[]*k1[] *XYZ[]];
+ Ers[] = AmplErs[] * Exp[I[]*k1r[]*XYZ[]];
+ Ets[] = AmplEts[] * Exp[I[]*k2[] *XYZ[]];
+ His[] = AmplHis[] * Exp[I[]*k1[] *XYZ[]];
+ Hrs[] = AmplHrs[] * Exp[I[]*k1r[]*XYZ[]];
+ Hts[] = AmplHts[] * Exp[I[]*k2[] *XYZ[]];
+ Eip[] = -1/(om0*ep0*epsr1[]) * k1[] /\ His[];
+ Erp[] = -1/(om0*ep0*epsr1[]) * k1r[] /\ Hrs[];
+ Etp[] = -1/(om0*ep0*epsr2[]) * k2[] /\ Hts[];
+
+ Ei[] = Ae*(Cos[psi0]*Eip[]-Sin[psi0]*Eis[]);
+ Er[] = Ae*(Cos[psi0]*Erp[]-Sin[psi0]*Ers[]);
+ Et[] = Ae*(Cos[psi0]*Etp[]-Sin[psi0]*Ets[]);
+ Hi[] = 1/(om0*mu0*mur[]) * k1[] /\ Ei[];
+ Hr[] = 1/(om0*mu0*mur[]) * k1r[] /\ Er[];
+ Ht[] = 1/(om0*mu0*mur[]) * k2[] /\ Et[];
+
+ E1[SurfIntTop] = Ei[]+Er[];
+ E1[Omega_super] = Ei[]+Er[];
+ E1[Omega_subs] = Et[];
+ E1[SurfIntBot] = Et[];
+ E1d[SurfIntTop] = Er[];
+ E1d[Omega_super] = Er[];
+ E1d[Omega_subs] = Et[];
+ E1d[SurfIntBot] = Et[];
+
+ H1[Omega_super] = Hi[]+Hr[];
+ H1[Omega_subs] = Ht[];
+ H1d[Omega_super] = Hr[];
+ H1d[Omega_subs] = Ht[];
+
+ source_vol_scat[] = (om0/cel)^2*(epsr[]-epsr_annex[])*E1[];
+ source_surf_tot[] = -2*I[]*om0*mu0*Vector[0,0,1] /\ Hi[];
+ // 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];
+ dephY[] = Exp[I[]*sky1[]*period_y];
+
+ // Fourier coefficients variables
+ Nb_ordre = 2*Nmax+1;
+ For i In {0:Nb_ordre-1}
+ For j In {0:Nb_ordre-1}
+ 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
- EndFor
- For i In {0:Nb_ordre-1}
- For j In {0:Nb_ordre-1}
- 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];
+ For i In {0:Nb_ordre-1}
+ For j In {0:Nb_ordre-1}
+ 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
- EndFor
-
-}
-
-Constraint {
- { Name Dirichlet; Type Assign;
- Case {
- { Region SurfDirichlet; Value 0.; }
- }
+
}
- { Name BlochX;
- Case {
- { Region SurfBlochXp; Type LinkCplx ; RegionRef SurfBlochXm;
- Coefficient dephX[]; Function Vector[$X-period_x,$Y,$Z] ; }
+
+ Constraint {
+ { Name Dirichlet; Type Assign;
+ Case {
+ { Region SurfDirichlet; Value 0.; }
+ }
}
- }
- { Name BlochY;
- Case {
- { Region SurfBlochYp; Type LinkCplx ; RegionRef SurfBlochYm;
- Coefficient dephY[]; Function Vector[$X-dys,$Y-dyc,$Z] ; }
+ { Name BlochX;
+ Case {
+ { Region SurfBlochXp; Type LinkCplx ; RegionRef SurfBlochXm;
+ Coefficient dephX[]; Function Vector[$X-period_x,$Y,$Z] ; }
+ }
}
- }
-}
-
-Jacobian {
- { Name JVol ; Case {{ Region All ; Jacobian Vol ; }}}
- { Name JSur ; Case {{ Region All ; Jacobian Sur ; }}}
- { Name JLin ; Case {{ Region All ; Jacobian Lin ; }}}
-}
-
-Integration {
- { Name I1 ;
- Case {
- { Type Gauss ;
- Case {
- { GeoElement Point ; NumberOfPoints 1 ; }
- { GeoElement Line ; NumberOfPoints 4 ; }
- { GeoElement Triangle ; NumberOfPoints 12 ; }
- { GeoElement Triangle2 ; NumberOfPoints 12 ; }
- { GeoElement Tetrahedron ; NumberOfPoints 16 ; }
- { GeoElement Tetrahedron2; NumberOfPoints 16 ; }
- }
+ { Name BlochY;
+ Case {
+ { Region SurfBlochYp; Type LinkCplx ; RegionRef SurfBlochYm;
+ Coefficient dephY[]; Function Vector[$X-dys,$Y-dyc,$Z] ; }
}
}
}
-}
-
-FunctionSpace {
- { Name Hcurl; Type Form1;
- BasisFunction {
- { Name sn; NameOfCoef un; Function BF_Edge; Support Region[{Omega}]; Entity EdgesOf[All]; }
- { Name sn2; NameOfCoef un2; Function BF_Edge_2E;Support Region[{Omega}]; Entity EdgesOf[All]; }
- If(oi==2)
- { Name sn3; NameOfCoef un3; Function BF_Edge_3F_b; Support Region[Omega]; Entity FacetsOf[Omega, Not SurfExcludeFacets]; }
- { Name sn4; NameOfCoef un4; Function BF_Edge_3F_c; Support Region[Omega]; Entity FacetsOf[Omega, Not SurfExcludeFacets]; }
- { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Omega]; Entity EdgesOf[Omega, Not SurfExcludeFacets]; }
- EndIf
- }
- Constraint {
- { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint BlochX; }
- { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint BlochY; }
- { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
- { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
- { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochY; }
- { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
- If (FlagLinkFacets==1)
- { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
- { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochY; }
- { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
- { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochY; }
- { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochX; }
- { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochY; }
- EndIf
- If(oi==2)
- { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
- { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
- { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
- EndIf
- }
+
+ Jacobian {
+ { Name JVol ; Case {{ Region All ; Jacobian Vol ; }}}
+ { Name JSur ; Case {{ Region All ; Jacobian Sur ; }}}
+ { Name JLin ; Case {{ Region All ; Jacobian Lin ; }}}
}
- { 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];}
+
+ Integration {
+ { Name I1 ;
+ Case {
+ { Type Gauss ;
+ Case {
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line ; NumberOfPoints 4 ; }
+ { GeoElement Triangle ; NumberOfPoints 12 ; }
+ { GeoElement Triangle2 ; NumberOfPoints 12 ; }
+ { GeoElement Tetrahedron ; NumberOfPoints 16 ; }
+ { GeoElement Tetrahedron2; NumberOfPoints 16 ; }
+ }
+ }
+ }
}
}
-}
-
-Formulation {
- { Name helmholtz_vector; Type FemEquation;
- Quantity {
- { Name u ; Type Local; NameOfSpace Hcurl; }
- { Name uz ; Type Local; NameOfSpace L2_lambda; }
+
+ FunctionSpace {
+ { Name Hcurl; Type Form1;
+ BasisFunction {
+ { Name sn; NameOfCoef un; Function BF_Edge; Support Region[{Omega,Gama}]; Entity EdgesOf[All]; }
+ { Name sn2; NameOfCoef un2; Function BF_Edge_2E;Support Region[{Omega,Gama}]; Entity EdgesOf[All]; }
+ If(oi==2)
+ { Name sn3; NameOfCoef un3; Function BF_Edge_3F_b; Support Region[{Omega,Gama}]; Entity FacetsOf[Omega, Not SurfExcludeFacets]; }
+ { Name sn4; NameOfCoef un4; Function BF_Edge_3F_c; Support Region[{Omega,Gama}]; Entity FacetsOf[Omega, Not SurfExcludeFacets]; }
+ { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[{Omega,Gama}]; Entity EdgesOf[Omega, Not SurfExcludeFacets]; }
+ EndIf
+ }
+ Constraint {
+ { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint BlochY; }
+ { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochY; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ If (FlagLinkFacets==1)
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochY; }
+ { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochY; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochY; }
+ EndIf
+ If(oi==2)
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ EndIf
+ }
}
- Equation{
- Galerkin { [-1/mur[] * Dof{Curl u} , {Curl u}]; In Omega ; Jacobian JVol; Integration I1; }
- Galerkin { [k0^2*epsr[] * Dof{u} , {u}]; In Omega ; Jacobian JVol; Integration I1; }
- If (FLAG_TOTAL==0)
- Galerkin { [source_vol_scat[] , {u}]; In Omega_source; Jacobian JVol; Integration I1; }
- 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;}
+ { 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];}
+ }
}
}
-}
-
-Resolution {
- { Name helmholtz_vector;
- System {
- { Name M; NameOfFormulation helmholtz_vector; Type ComplexValue; }
+
+ Formulation {
+ { Name helmholtz_vector; Type FemEquation;
+ Quantity {
+ { 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; }
+ Galerkin { [k0^2*epsr[] * Dof{u} , {u}]; In Omega ; Jacobian JVol; Integration I1; }
+ If (FLAG_TOTAL==0)
+ Galerkin { [source_vol_scat[] , {u}]; In Omega_source; Jacobian JVol; Integration I1; }
+ 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;}
+ }
}
- Operation {
- CreateDir[Str[myDir]];
- Generate[M];
- Solve[M]; //SaveSolution[M];
+ }
+
+ Resolution {
+ { Name helmholtz_vector;
+ System {
+ { Name M; NameOfFormulation helmholtz_vector; Type ComplexValue; }
+ }
+ Operation {
+ CreateDir[Str[myDir]];
+ Generate[M];
+ Solve[M]; //SaveSolution[M];
+ }
}
}
-}
-
-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; } } }
- { Name E1 ; Value { Local { [ E1[] ]; In Omega; Jacobian JVol; } } }
- { Name lambda_step ; Value { Local { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
-
- If (FLAG_TOTAL==0)
- { 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 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; } } }
- For k In {2:6}
- { 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*dyc) ] ; In Scat ; Integration I1 ; Jacobian JVol ; } } }
- Else
- { Name Etot ; Value { Local { [ {u} ]; In Omega; Jacobian JVol; } } }
- { Name Htot ; Value { Local { [ -I[]/(mur[]*mu0*om0)*{Curl u}]; In Omega; Jacobian JVol; } } }
- { Name Poy_tot; Value { Local { [ 0.5*Re[Cross[{u} , Conj[ -I[]/(mur[]*mu0*om0)*{Curl u}]]] ]; In Omega; Jacobian JVol; } } }
- { Name Poy_ref; Value { Local { [ 0.5*Re[Cross[{u}-Ei[], Conj[-Hi[]-I[]/(mur[]*mu0*om0)*{Curl u}]]] ]; In Omega; Jacobian JVol; } } }
- For k In {2:6}
- { Name Abs_L~{k} ; Value { Integral { [ ep0*om0 * 0.5*Im[CompXX[epsr[]]]*(SquNorm[{u}]) / (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}]) / (Pinc*period_x*dyc) ] ; In Scat ; Integration I1 ; Jacobian JVol ; } } }
- EndIf
-
- 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_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 ; } } }
+
+ PostProcessing {
+ { Name postpro_helmholtz_vector; NameOfFormulation helmholtz_vector; NameOfSystem M;
+ Quantity {
+ { Name u ; Value { Local { [ {u} ]; In Omega; Jacobian JVol; } } }
+ { Name CompZu ; Value { Local { [ CompZ[{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; } } }
+ { Name E1 ; Value { Local { [ E1[] ]; In Omega; Jacobian JVol; } } }
+ { Name lambda_step ; Value { Local { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
+
+ If (FLAG_TOTAL==0)
+ { 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 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; } } }
+ For k In {2:6}
+ { 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*dyc) ] ; In Scat ; Integration I1 ; Jacobian JVol ; } } }
+ Else
+ { Name Etot ; Value { Local { [ {u} ]; In Omega; Jacobian JVol; } } }
+ { Name Htot ; Value { Local { [ -I[]/(mur[]*mu0*om0)*{Curl u}]; In Omega; Jacobian JVol; } } }
+ { Name Poy_tot; Value { Local { [ 0.5*Re[Cross[{u} , Conj[ -I[]/(mur[]*mu0*om0)*{Curl u}]]] ]; In Omega; Jacobian JVol; } } }
+ { Name Poy_ref; Value { Local { [ 0.5*Re[Cross[{u}-Ei[], Conj[-Hi[]-I[]/(mur[]*mu0*om0)*{Curl u}]]] ]; In Omega; Jacobian JVol; } } }
+ For k In {2:6}
+ { Name Abs_L~{k} ; Value { Integral { [ ep0*om0 * 0.5*Im[CompXX[epsr[]]]*(SquNorm[{u}]) / (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}]) / (Pinc*period_x*dyc) ] ; In Scat ; Integration I1 ; Jacobian JVol ; } } }
+ EndIf
+
+ 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_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_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_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 eff_t2~{i}~{j} ; Value { Term{ Type Global; [
- gammat~{i}~{j}[]/(-k1z[]*Cos[xsi]^2) * ( SquNorm[$int_x_t~{i}~{j}]+
- SquNorm[$int_y_t~{i}~{j}]+
- SquNorm[$int_z_t~{i}~{j}] ) ] ; In SurfIntBot ; } } }
- { Name eff_r2~{i}~{j} ; Value { Term{ Type Global; [
- gammar~{i}~{j}[]/(-k1z[]*Cos[xsi]^2) * ( SquNorm[$int_x_r~{i}~{j}]+
- SquNorm[$int_y_r~{i}~{j}]+
- 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 ; } } }
+ { 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_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_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 eff_t2~{i}~{j} ; Value { Term{ Type Global; [
+ gammat~{i}~{j}[]/(-k1z[]*Cos[xsi]^2) * ( SquNorm[$int_x_t~{i}~{j}]+
+ SquNorm[$int_y_t~{i}~{j}]+
+ SquNorm[$int_z_t~{i}~{j}] ) ] ; In SurfIntBot ; } } }
+ { Name eff_r2~{i}~{j} ; Value { Term{ Type Global; [
+ gammar~{i}~{j}[]/(-k1z[]*Cos[xsi]^2) * ( SquNorm[$int_x_r~{i}~{j}]+
+ SquNorm[$int_y_r~{i}~{j}]+
+ 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
- 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 ; } } }
- { Name et_specular ; Value { Term{Type Global; [
- Vector[$int_x_t~{ispecular}~{jspecular},
- $int_y_t~{ispecular}~{jspecular},
- $int_z_t~{ispecular}~{jspecular}] ] ; In SurfIntBot ; } } }
-
- // Project er_specular on the (s,p) basis
- { Name rp ; Value { Term{ Type Global; [ ppol_r[] * $er_specular] ; In SurfIntTop ; } } }
- { Name rs ; Value { Term{ Type Global; [ spol[] * $er_specular] ; In SurfIntTop ; } } }
- { Name tp ; Value { Term{ Type Global; [ ppol_t[] * $et_specular] ; In SurfIntBot ; } } }
- { Name ts ; Value { Term{ Type Global; [ spol[] * $et_specular] ; In SurfIntBot ; } } }
-
-}
+ // 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 ; } } }
+ { Name et_specular ; Value { Term{Type Global; [
+ Vector[$int_x_t~{ispecular}~{jspecular},
+ $int_y_t~{ispecular}~{jspecular},
+ $int_z_t~{ispecular}~{jspecular}] ] ; In SurfIntBot ; } } }
+
+ // Project er_specular on the (s,p) basis
+ { Name rp ; Value { Term{ Type Global; [ ppol_r[] * $er_specular] ; In SurfIntTop ; } } }
+ { Name rs ; Value { Term{ Type Global; [ spol[] * $er_specular] ; In SurfIntTop ; } } }
+ { Name tp ; Value { Term{ Type Global; [ ppol_t[] * $et_specular] ; In SurfIntBot ; } } }
+ { Name ts ; Value { Term{ Type Global; [ spol[] * $et_specular] ; In SurfIntBot ; } } }
+
}
-}
-
-PostOperation {
- { Name postop_helmholtz_vector; NameOfPostProcessing postpro_helmholtz_vector ;
- Operation {
- 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"];
- Else
- Print [ u , OnElementsOf Omega, File StrCat[myDir,"Etot.pos"]];
+ }
+ }
+
+ PostOperation {
+ { Name postop_helmholtz_vector; NameOfPostProcessing postpro_helmholtz_vector ;
+ Operation {
+
+ Print [ uz , OnElementsOf SurfIntTop, File StrCat[myDir,"uz_ZP.pos"]];
+ Print [ uz , OnElementsOf SurfIntBot, File StrCat[myDir,"uz_ZM.pos"]];
+
+ Print [ epsr_xx , OnElementsOf Omega, File StrCat[myDir,"epsr_xx.pos"]];
+ Print [ u , OnElementsOf Omega, File StrCat[myDir,"Edif.pos"]];
+ Print [ Etot , OnElementsOf Omega, File StrCat[myDir,"Etot.pos"]];
+ Print [ CompZu , OnPlane { {-period_x/2,-period_y/2,hh_L_6+0.5*nm} { period_x/2,-period_y/2,hh_L_6+0.5*nm} {-period_x/2, period_y/2,hh_L_6+0.5*nm} } {npts_interpX,npts_interpY} , File StrCat[myDir,"u_cut_ZM.pos"]];
+ Print [ CompZu , OnPlane { {-period_x/2,-period_y/2,hh_L_1+thick_L_1-0.5*nm} { period_x/2,-period_y/2,hh_L_1+thick_L_1-0.5*nm} {-period_x/2, period_y/2,hh_L_1+thick_L_1-0.5*nm} } {npts_interpX,npts_interpY} , File StrCat[myDir,"u_cut_ZP.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"];
+ Else
+ Print [ u , OnElementsOf Omega, File StrCat[myDir,"Etot.pos"]];
+ EndIf
EndIf
- EndIf
- If (FlagOutEtotFull==1)
- If (Flag_interp_cubic==1)
- Print [ Etot , OnBox { {-period_x/2,-period_y/2,hh_L_6} {period_x/2,-period_y/2,hh_L_6} {-period_x/2,period_y/2,hh_L_6} {-period_x/2,-period_y/2,hh_L_1+thick_L_1} } {npts_interpX,npts_interpY,npts_interpZTot} , File StrCat[myDir,"Etot_grid.pos"], Name "Etot_grid"];
- Else
- Print [ Etot , OnElementsOf Omega_plot, File StrCat[myDir,"Etot.pos"]];
+ If (FlagOutEtotFull==1)
+ If (Flag_interp_cubic==1)
+ Print [ Etot , OnBox { {-period_x/2,-period_y/2,hh_L_6} {period_x/2,-period_y/2,hh_L_6} {-period_x/2,period_y/2,hh_L_6} {-period_x/2,-period_y/2,hh_L_1+thick_L_1} } {npts_interpX,npts_interpY,npts_interpZTot} , File StrCat[myDir,"Etot_grid.pos"], Name "Etot_grid"];
+ Else
+ Print [ Etot , OnElementsOf Omega_plot, File StrCat[myDir,"Etot.pos"]];
+ EndIf
EndIf
- EndIf
- If (FlagOutPoyFull==1)
- If (Flag_interp_cubic==1)
- Print [ Poy_tot , OnBox { {-period_x/2,-period_y/2,hh_L_6} {period_x/2,-period_y/2,hh_L_6} {-period_x/2,period_y/2,hh_L_6} {-period_x/2,-period_y/2,hh_L_1+thick_L_1} } {npts_interpX,npts_interpY,npts_interpZTot} , File StrCat[myDir,"Poy_tot_grid.pos"], Name "Poy_tot_grid"];
- Else
- Print [ Poy_tot , OnElementsOf Omega_plot, File StrCat[myDir,"Poy_tot.pos"]];
+ If (FlagOutPoyFull==1)
+ If (Flag_interp_cubic==1)
+ Print [ Poy_tot , OnBox { {-period_x/2,-period_y/2,hh_L_6} {period_x/2,-period_y/2,hh_L_6} {-period_x/2,period_y/2,hh_L_6} {-period_x/2,-period_y/2,hh_L_1+thick_L_1} } {npts_interpX,npts_interpY,npts_interpZTot} , File StrCat[myDir,"Poy_tot_grid.pos"], Name "Poy_tot_grid"];
+ Else
+ Print [ Poy_tot , OnElementsOf Omega_plot, File StrCat[myDir,"Poy_tot.pos"]];
+ EndIf
EndIf
- EndIf
- If (FlagOutEscaCuts==1)
- Print [ u , OnPlane { {-period_x/2,0,hh_L_6-PML_bot} {period_x/2,0,hh_L_6-PML_bot} {-period_x/2,0,hh_L_1+thick_L_1+PML_top} } {npts_interpX,npts_interpZSca} , File StrCat[myDir,"u_cut_Y=0.pos"], Name "u_cut_Y=0"];
- Print [ u , OnPlane { {0,-period_y/2,hh_L_6-PML_bot} {0,period_y/2,hh_L_6-PML_bot} {0,-period_y/2,hh_L_1+thick_L_1+PML_top} } {npts_interpY,npts_interpZSca} , File StrCat[myDir,"u_cut_X=0.pos"], Name "u_cut_X=0"];
- EndIf
- If (FlagOutEtotCuts==1)
- Print [ Etot , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZTot} , File StrCat[myDir,"Etot_cut_Y=0.pos"], Name "Etot_cut_Y=0"];
- Print [ Etot , OnPlane { {0,-period_y/2,hh_L_6} {0,period_y/2,hh_L_6} {0,-period_y/2,hh_L_1+thick_L_1} } {npts_interpY,npts_interpZTot} , File StrCat[myDir,"Etot_cut_X=0.pos"], Name "Etot_cut_X=0"];
- EndIf
- If (FlagOutHtotCuts==1)
- Print [ Htot , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZTot} , File StrCat[myDir,"Htot_cut_Y=0.pos"], Name "Htot_cut_Y=0"];
- Print [ Htot , OnPlane { {0,-period_y/2,hh_L_6} {0,period_y/2,hh_L_6} {0,-period_y/2,hh_L_1+thick_L_1} } {npts_interpY,npts_interpZTot} , File StrCat[myDir,"Htot_cut_X=0.pos"], Name "Htot_cut_X=0"];
- EndIf
- If (FlagOutPoyCut==1)
- Print [ Poy_tot , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZTot} , File StrCat[myDir,"Poy_tot_cut_Y=0.pos"], Name "Poy_tot_cut_Y=0"];
- Print [ Poy_tot , OnPlane { {0,-period_y/2,hh_L_6} {0,period_y/2,hh_L_6} {0,-period_y/2,hh_L_1+thick_L_1} } {npts_interpY,npts_interpZTot} , File StrCat[myDir,"Poy_tot_cut_X=0.pos"], Name "Poy_tot_cut_X=0"];
- EndIf
-
- 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_checkpoyX-1,npts_checkpoyY-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_checkpoyX-1,npts_checkpoyY-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_checkpoyX-1,npts_checkpoyY-1} , File StrCat[myDir,"Poy_inc_gd.pos"], Format Table];
-
- 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
- Print[ Abs_scat[Scat] , OnGlobal, File > StrCat[myDir,"temp-Q_scat.txt"], Format Table ];
-
- For i In {0:Nb_ordre-1}
- 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];
+ If (FlagOutEscaCuts==1)
+ Print [ u , OnPlane { {-period_x/2,0,hh_L_6-PML_bot} {period_x/2,0,hh_L_6-PML_bot} {-period_x/2,0,hh_L_1+thick_L_1+PML_top} } {npts_interpX,npts_interpZSca} , File StrCat[myDir,"u_cut_Y=0.pos"], Name "u_cut_Y=0"];
+ Print [ u , OnPlane { {0,-period_y/2,hh_L_6-PML_bot} {0,period_y/2,hh_L_6-PML_bot} {0,-period_y/2,hh_L_1+thick_L_1+PML_top} } {npts_interpY,npts_interpZSca} , File StrCat[myDir,"u_cut_X=0.pos"], Name "u_cut_X=0"];
+ EndIf
+ If (FlagOutEtotCuts==1)
+ Print [ Etot , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZTot} , File StrCat[myDir,"Etot_cut_Y=0.pos"], Name "Etot_cut_Y=0"];
+ Print [ Etot , OnPlane { {0,-period_y/2,hh_L_6} {0,period_y/2,hh_L_6} {0,-period_y/2,hh_L_1+thick_L_1} } {npts_interpY,npts_interpZTot} , File StrCat[myDir,"Etot_cut_X=0.pos"], Name "Etot_cut_X=0"];
+ Print [ E1 , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZTot} , File StrCat[myDir,"E1_cut_Y=0.pos"], Name "E1_cut_Y=0"];
+ Print [ E1 , OnPlane { {0,-period_y/2,hh_L_6} {0,period_y/2,hh_L_6} {0,-period_y/2,hh_L_1+thick_L_1} } {npts_interpY,npts_interpZTot} , File StrCat[myDir,"E1_cut_X=0.pos"], Name "E1_cut_X=0"];
+ EndIf
+ If (FlagOutHtotCuts==1)
+ Print [ Htot , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZTot} , File StrCat[myDir,"Htot_cut_Y=0.pos"], Name "Htot_cut_Y=0"];
+ Print [ Htot , OnPlane { {0,-period_y/2,hh_L_6} {0,period_y/2,hh_L_6} {0,-period_y/2,hh_L_1+thick_L_1} } {npts_interpY,npts_interpZTot} , File StrCat[myDir,"Htot_cut_X=0.pos"], Name "Htot_cut_X=0"];
+ EndIf
+ If (FlagOutPoyCut==1)
+ Print [ Poy_tot , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZTot} , File StrCat[myDir,"Poy_tot_cut_Y=0.pos"], Name "Poy_tot_cut_Y=0"];
+ Print [ Poy_tot , OnPlane { {0,-period_y/2,hh_L_6} {0,period_y/2,hh_L_6} {0,-period_y/2,hh_L_1+thick_L_1} } {npts_interpY,npts_interpZTot} , File StrCat[myDir,"Poy_tot_cut_X=0.pos"], Name "Poy_tot_cut_X=0"];
+ EndIf
+
+ 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_checkpoyX-1,npts_checkpoyY-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_checkpoyX-1,npts_checkpoyY-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_checkpoyX-1,npts_checkpoyY-1} , File StrCat[myDir,"Poy_inc_gd.pos"], Format Table];
+
+ 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
- EndFor
-
- For i In {0:Nb_ordre-1}
- For j In {0:Nb_ordre-1}
- 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 ];
+ Print[ Abs_scat[Scat] , OnGlobal, File > StrCat[myDir,"temp-Q_scat.txt"], Format Table ];
+
+ For i In {0:Nb_ordre-1}
+ 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
- EndFor
- Print[ er_specular[SurfIntTop] , OnRegion SurfIntTop, StoreInVariable $er_specular, Format Table];
- Print[ et_specular[SurfIntBot] , OnRegion SurfIntBot, StoreInVariable $et_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 ];
- Print[ tp[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir,"tp.txt"], Format Table ];
- Print[ ts[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir,"ts.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 ];
+
+ For i In {0:Nb_ordre-1}
+ For j In {0:Nb_ordre-1}
+ 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
- EndFor
- Print[ lambda_step, OnPoint{0,0,0}, Format Table, File > StrCat[myDir, "temp_lambda_step.txt"], SendToServer "GetDP/Lambda_step" ] ;
+ Print[ er_specular[SurfIntTop] , OnRegion SurfIntTop, StoreInVariable $er_specular, Format Table];
+ Print[ et_specular[SurfIntBot] , OnRegion SurfIntBot, StoreInVariable $et_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 ];
+ Print[ tp[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir,"tp.txt"], Format Table ];
+ Print[ ts[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir,"ts.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 ];
+ EndFor
+ EndFor
+ Print[ lambda_step, OnPoint{0,0,0}, Format Table, File > StrCat[myDir, "temp_lambda_step.txt"], SendToServer "GetDP/Lambda_step" ] ;
+ }
}
}
-}
-
-DefineConstant[
- R_ = {"helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
- 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}
-];
+
+ DefineConstant[
+ R_ = {"helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+ 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}
+ ];
+
\ No newline at end of file
diff --git a/DiffractionGratings/grating3D_data_half_ellipsoid.geo b/DiffractionGratings/grating3D_data_half_ellipsoid.geo
index 07906f9a4a7e90110c3362b25b98dbc1f7e359ef..38f5932acb3551f1998d73cd8eaf1931259e00d7 100644
--- a/DiffractionGratings/grating3D_data_half_ellipsoid.geo
+++ b/DiffractionGratings/grating3D_data_half_ellipsoid.geo
@@ -1,4 +1,4 @@
-nm = 1000;
+nm = 1;
pp1 = "1Incident Plane Wave";
pp2 = "2Layers Thicknesses";
pp3 = "3Scatterer Properties";
@@ -11,8 +11,8 @@ DefineConstant[
thetadeg = {40 , Name StrCat[pp1,"/2theta0 [deg]"]},
phideg = {36 , Name StrCat[pp1,"/3phi0 [deg]"]},
psideg = {72 , Name StrCat[pp1,"/4psi0 [deg]"]},
- period_x = {250 , Name StrCat[pp2,"/1X period [nm]"]},
- period_y = {250 , Name StrCat[pp2,"/2Y period [nm]"]},
+ period_x = {473 , Name StrCat[pp2,"/1X period [nm]"]},
+ period_y = {322 , Name StrCat[pp2,"/2Y period [nm]"]},
thick_L_1 = {50 , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
thick_L_2 = {50 , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
thick_L_3 = {100 , Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
@@ -21,29 +21,29 @@ DefineConstant[
thick_L_6 = {50 , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
xsideg = {0 , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
- tag_geom = { 4 , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
- rx = {107 , Name StrCat[pp3,"/1rx"]},
- ry = {47 , Name StrCat[pp3,"/2ry"]},
- rz = {40 , Name StrCat[pp3,"/3rz"]},
- flag_mat_scat = { 0 , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
- eps_re_Scat = {-2.23, Name StrCat[pp3,"/7eps_re_Scat"]},
- eps_im_Scat = { 3.89, Name StrCat[pp3,"/8eps_im_Scat"]},
+ tag_geom = { 4 , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
+ rx = {107 , Name StrCat[pp3,"/1rx"]},
+ ry = {47 , Name StrCat[pp3,"/2ry"]},
+ rz = {40 , Name StrCat[pp3,"/3rz"]},
+ flag_mat_scat = { 4 , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ eps_re_Scat = { 9 , Name StrCat[pp3,"/7eps_re_Scat"]},
+ eps_im_Scat = { 1 , Name StrCat[pp3,"/8eps_im_Scat"]},
flag_mat_1 = { 0 , Name StrCat[pp4,"/1Layer 1"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
flag_mat_2 = { 0 , Name StrCat[pp4,"/2Layer 2"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
flag_mat_3 = { 0 , Name StrCat[pp4,"/3Layer 3"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
- flag_mat_4 = { 0 , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
- flag_mat_5 = { 0 , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
- flag_mat_6 = { 0 , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_4 = { 1 , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_5 = { 1 , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_6 = { 1 , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
eps_re_L_1 = {1 , Name StrCat[pp4,"/layer 1: real part of relative permittivity"]},
eps_im_L_1 = {0 , Name StrCat[pp4,"/layer 1: imag part of relative permittivity"]},
eps_re_L_2 = {1 , Name StrCat[pp4,"/layer 2: real part of relative permittivity"]},
eps_im_L_2 = {0 , Name StrCat[pp4,"/layer 2: imag part of relative permittivity"]},
eps_re_L_3 = {1 , Name StrCat[pp4,"/layer 3: real part of relative permittivity"]},
eps_im_L_3 = {0 , Name StrCat[pp4,"/layer 3: imag part of relative permittivity"]},
- eps_re_L_4 = {1 , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
+ eps_re_L_4 = {4 , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
eps_im_L_4 = {0 , Name StrCat[pp4,"/layer 4: imag part of relative permittivity"]},
- eps_re_L_5 = {1 , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
+ eps_re_L_5 = {4 , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
eps_im_L_5 = {0 , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
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"]},
@@ -51,9 +51,9 @@ DefineConstant[
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]"]},
- PML_top = {lambda0 , Name StrCat[pp5,"/4PML top thickness [nm]"]},
- PML_bot = {lambda0 , Name StrCat[pp5,"/5PML bot thickness [nm]"]},
+ lc_scat = {1 , Name StrCat[pp5,"/2metal mesh size [nm]"]},
+ PML_top = {800 , Name StrCat[pp5,"/4PML top thickness [nm]"]},
+ PML_bot = {800 , Name StrCat[pp5,"/5PML bot thickness [nm]"]},
Nmax = {1 , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
refine_mesh_L_1= {1 , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
refine_mesh_L_2= {1 , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
@@ -64,11 +64,11 @@ DefineConstant[
FlagLinkFacets = {0 , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
PML_TYPE = {0 , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
- InterpSampling = { 10 , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
+ InterpSampling = { 2 , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
Flag_interp_cubic = { 0 , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
FlagOutEtotCuts = { 1 , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
FlagOutHtotCuts = { 0 , Name StrCat[pp6,"/3Output Total Magnetic Field cuts?"] , Choices {0,1} },
- FlagOutEscaCuts = { 0 , Name StrCat[pp6,"/4Output Scattered Electric Field cuts?"] , Choices {0,1} },
+ FlagOutEscaCuts = { 1 , Name StrCat[pp6,"/4Output Scattered Electric Field cuts?"] , Choices {0,1} },
FlagOutPoyCut = { 0 , 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} },
diff --git a/DiffractionGratings/grating3D_parallel_Mmatrix.sh b/DiffractionGratings/grating3D_parallel_Mmatrix.sh
index 887d864165e8fa66fed6838ae31e15c1725e05ba..123ccc32d66783111a005276edc3df152a2c414c 100644
--- a/DiffractionGratings/grating3D_parallel_Mmatrix.sh
+++ b/DiffractionGratings/grating3D_parallel_Mmatrix.sh
@@ -4,10 +4,10 @@ export OMP_NUM_THREADS=1
export NPROC=96
# export nb_lam=1
# export nb_phi=1
-export nb_lam=100
-export nb_phi=59
+export nb_lam=20
+export nb_phi=20
export lambda_min=400
-export lambda_max=1600
+export lambda_max=1200
export FLAG_TOTAL=0
export GRATING_CASE="half_ellipsoid"
export myDir="res_Matrix_nb_lam"$nb_lam"_nb_phi"$nb_phi"_total"$FLAG_TOTAL
@@ -20,7 +20,12 @@ myfunc() {
local mysubDir=$myDir/run_lam$1_phi$2_psi$3
mkdir $mysubDir
cp grating3D.msh grating3D.pro grating3D_data_$GRATING_CASE.geo grating3D_data.geo grating3D_materials.pro $mysubDir
- local lam=$(echo "scale=10;400+400/($nb_lam-1)*$1" | bc )
+ if (( $nb_lam == 1 ))
+ then
+ local lam=$(echo "scale=10;$lambda_min" | bc )
+ else
+ local lam=$(echo "scale=10;$lambda_min+($lambda_max-$lambda_min)/($nb_lam-1)*$1" | bc )
+ fi
local phi=$(echo "scale=10;360/($nb_phi)*$2" | bc )
local psi=$(echo "scale=10;90*$3" | bc )
cd $mysubDir
diff --git a/DiffractionGratings/grating3D_postplot_Mmatrix.py b/DiffractionGratings/grating3D_postplot_Mmatrix.py
index d3a0975feec718638789c90e8893e50e80e73ebd..4b7269304160f58c377882e8391678d54e8738af 100644
--- a/DiffractionGratings/grating3D_postplot_Mmatrix.py
+++ b/DiffractionGratings/grating3D_postplot_Mmatrix.py
@@ -1,6 +1,7 @@
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
+import scipy.constants as scc
plt.rcParams.update({"text.usetex": True, "font.family": "serif"})
def add_colorbar(mappable):
@@ -22,10 +23,10 @@ def load_getdp_integral(fname):
else:
return temp[:,1]+1j*temp[:,2]
-nb_lam = 200
-nb_phi = 60
+nb_lam = 20
+nb_phi = 21
FLAG_TOT = 0
-respath = 'res_Matrix_nb_lam200_nb_phi59_total0/'
+respath = 'res_Matrix_nb_lam%g_nb_phi%g_total0/'%(nb_lam,nb_phi-1)
rpp = np.zeros((nb_lam,nb_phi),dtype=complex)
rps = np.zeros((nb_lam,nb_phi),dtype=complex)
rsp = np.zeros((nb_lam,nb_phi),dtype=complex)
@@ -42,8 +43,9 @@ effr2_sin = np.zeros((9,nb_lam,nb_phi),dtype=complex)
Qscat_sin = np.zeros((nb_lam,nb_phi))
-tab_lam = np.linspace(350,1600,nb_lam)
+tab_lam = np.linspace(495,1200,nb_lam)
tab_phi = np.linspace(0,360,nb_phi)
+tab_hnu = 2*scc.pi*scc.c/(tab_lam*1e-9/scc.hbar*scc.eV)
M = np.zeros((4,4,len(tab_lam),len(tab_phi)),dtype=complex)
@@ -89,6 +91,10 @@ T1_pin = np.sum(efft1_pin,axis=0)
T2_pin = np.sum(efft2_pin,axis=0)
R1_pin = np.sum(effr1_pin,axis=0)
R2_pin = np.sum(effr2_pin,axis=0)
+T00_sin = efft2_sin[4,:,:]
+R00_sin = effr2_sin[4,:,:]
+T00_pin = efft2_pin[4,:,:]
+R00_pin = effr2_pin[4,:,:]
TOT1_pin = T1_pin+R1_pin+Qscat_pin
TOT2_pin = T2_pin+R2_pin+Qscat_pin
TOT1_sin = T1_sin+R1_sin+Qscat_sin
@@ -137,22 +143,35 @@ for form in range(2):
axes[2,1].title.set_text('effr1 pin')
plt.savefig('BALANCE_FLAG_TOT%g.jpg'%FLAG_TOT)
+fig, axes = plt.subplots(3, 2, figsize=(12,12))
+L,P = np.meshgrid(tab_lam,tab_phi,indexing='ij')
+ax=axes[0,0];zplot = ax.pcolormesh(L,P,T00_sin.real );add_colorbar(zplot);ax.title.set_text('T00 sin')
+ax=axes[1,0];zplot = ax.pcolormesh(L,P,R00_sin.real );add_colorbar(zplot);ax.title.set_text('R00 sin')
+ax=axes[2,0];zplot = ax.pcolormesh(L,P,Qscat_sin.real);add_colorbar(zplot);ax.title.set_text('Abs sin')
+ax=axes[0,1];zplot = ax.pcolormesh(L,P,T00_pin.real );add_colorbar(zplot);ax.title.set_text('T00 pin')
+ax=axes[1,1];zplot = ax.pcolormesh(L,P,R00_pin.real );add_colorbar(zplot);ax.title.set_text('R00 pin')
+ax=axes[2,1];zplot = ax.pcolormesh(L,P,Qscat_pin.real);add_colorbar(zplot);ax.title.set_text('Abs pin')
+plt.savefig('BALANCE_new_code.jpg')
+
+
+
fig, axes = plt.subplots(4, 4, subplot_kw=dict(projection='polar') ,figsize=(12,12))
-flag_lam = True
-rlabel = r"$\lambda_0$ (nm)"
-which_r = tab_lam
+flag_lam = False
+rlabel = r"$\hbar \nu$" if not flag_lam else r"$\lambda_0$"
+which_r = tab_hnu if not flag_lam else tab_lam
which_orig = 0 #if not flag_lam else 350
for i in range(4):
for j in range(4):
ax=axes[i,j]
if i!=0 or j!=0:
- sm=ax.contourf(tab_phi*np.pi/180,which_r,M[i,j]/M[0,0],cmap=plt.cm.bwr,levels=30)
+ sm=ax.contourf(tab_phi*np.pi/180,which_r,M[i,j]/M[0,0],cmap=plt.cm.bwr,levels=30,vmin=-1,vmax=1)
+ # sm=ax.contourf(tab_phi*np.pi/180,which_r,M[i,j]/M[0,0],cmap=plt.cm.bwr,levels=30)
ax.text(0. , 1. , r"$M_{%d%d}$"%(i+1,j+1) , fontsize=16 , transform=ax.transAxes)
ax.set_xticks([])
ax.set_yticks([])
- cbar = plt.colorbar(sm, ax=ax, fraction=0.046, pad=0.04)
- cbar.ax.locator_params(nbins=5)
- cbar.ax.tick_params(labelsize=14)
+ # cbar = plt.colorbar(sm, ax=ax, fraction=0.046, pad=0.04)
+ # cbar.ax.locator_params(nbins=5)
+ # cbar.ax.tick_params(labelsize=14)
else:
p00 = ax.contourf(tab_phi*np.pi/180,which_r,M[i,j]/M[0,0]-1,cmap=plt.cm.bwr,vmin=-1,vmax=1)
ax.xaxis.label.set_color('C0') #setting up X-axis label color to yellow