diff --git a/ECE_full_wave_E-H_monopole/geo/PhysicalRegions.geo b/ECE_full_wave_E-H_monopole/geo/PhysicalRegions.geo
new file mode 100644
index 0000000000000000000000000000000000000000..fdb9e591c09a93b66fe17923c019bdf2935d3699
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/geo/PhysicalRegions.geo
@@ -0,0 +1,27 @@
+//=================================================
+// Physical regions
+//=================================================
+
+// surfaces
+Physical Surface(MONOPOLE) = {surfMonopole[]};
+Physical Surface(AIR) = surfAir;
+
+// boundaries
+Physical Line(TERMINAL) = {feedMonople[]}; // Feeding
+
+Physical Line(CUT) = {lbox[{0}]};
+Physical Line(AXIS) = {lbox[{3}],axisMonopole[]};
+
+If (Flag_AnalysisTypeExtBnd == 1)
+ Physical Line(GROUND) = {lbox[{1}],lbox[{2}]};
+Else
+ Physical Line(SURFAIRINF) = interfaceAirPML;
+ Physical Line(GROUND) = {lbox[{1}]};
+EndIf
+
+
+// For aestetics
+Mesh.Light = 0; // so that to see mesh colors in the coax part
+
+Recursive Color SkyBlue { Surface{surfAir};}
+Recursive Color Black { Surface{surfMonopole};}
\ No newline at end of file
diff --git a/ECE_full_wave_E-H_monopole/geo/monopole2D_cap_pml.geo b/ECE_full_wave_E-H_monopole/geo/monopole2D_cap_pml.geo
new file mode 100644
index 0000000000000000000000000000000000000000..ad75f9bded36235d1cb37cb94c66e3f59f0be9c9
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/geo/monopole2D_cap_pml.geo
@@ -0,0 +1,147 @@
+Include "../monopole_data.pro";
+
+fac = s;
+flag = 1; // if 0, extruded mesh becomes free in the monopole
+_use_structured_pml = 0 ;
+
+lcdAntennaZero = a/4;
+lcdAntennaA = lcdAntennaZero;
+
+If (_use_fineMesh == 1)
+ lambdaMin = c0/fmax;
+ lcd = fac*lambdaMin/nopPerLambdaAir;
+ lcair = lcd;
+
+ nopL = flag*(Ceil[L/lcdAntennaZero]); // to make the antenna very fine along L
+Else
+ lcdAntennaZero = 4*lcdAntennaZero;
+ lambdaMin = c0/fmax;
+ lcair = fac*lambdaMin/nopPerLambdaAir;
+
+ delta = Sqrt(2.0*nuMonopole/(2*Pi*Freq*sigmaMonopole));
+ lcdAntennaA = fac*delta/nbDelta;
+
+ lcair1 = lambdaMin/10;
+ nopL = fac*6*flag*(Ceil[L/lcair1]);
+EndIf
+nopL = (nopL < 10) ? 10 : nopL;
+
+Printf("nopL %g",nopL);
+Printf("lcdAntennaA %g lcdAntennaZero %g ratio %g", lcdAntennaA, lcdAntennaZero, lcdAntennaZero/lcdAntennaA);
+
+
+//=================================================
+// Monopole
+//=================================================
+
+d = 0; // no gap
+yd = (d >= 0)? d:0.;
+p0 = newp ; Point(p0) = {0, yd, 0, lcdAntennaZero};
+p1 = newp ; Point(p1) = {a, yd, 0, lcdAntennaA};
+LineMonopoleBottom = newl ; Line(LineMonopoleBottom) = {p0,p1};
+
+surf[] = Extrude {0, L, 0} {
+ Line{LineMonopoleBottom}; Layers{nopL} ;
+};
+surfMonopole[0] = surf[1] ; LineMonopoleTop = surf[0] ;
+
+p_[] = Boundary{Line{LineMonopoleTop};};
+
+interfaceMonopoleAir[] = Boundary{ Surface{surfMonopole[]};} ;
+
+axisMonopole[] = interfaceMonopoleAir[{3}] ;
+feedMonople[] = interfaceMonopoleAir[{0}] ;
+interfaceMonopoleAir[] -= axisMonopole[];
+interfaceMonopoleAir[] -= feedMonople[];
+
+//=================================================
+// Air and PML
+//=================================================
+
+pZero = newp; Point(pZero) = { 0, 0, 0, lcdAntennaZero};
+
+pp[]+=p1;
+pp[]+=newp; Point(newp) = { b, 0, 0, lcdAntennaZero};
+
+// reduce air around
+cc = 1; // 4
+pp[]+=newp; Point(newp) = { xb/cc, 0, 0, lcair};
+pp[]+=newp; Point(newp) = { xb/cc, L, 0, lcair};
+pp[]+=newp; Point(newp) = { 0, xb/cc+L, 0, lcair};
+pp[]+=p_[0];
+
+For k In {0:2}
+ lbox[]+=newl; Line(newl) = {pp[k],pp[k+1]};
+EndFor
+lbox[]+=newl; Circle(newl) = {pp[3],p_[0],pp[4]};
+lbox[]+=newl; Line(newl) = {pp[4],pp[5]};
+interfaceAirPML[] = lbox[{2,3}]; // it will be the inf boundary if no PML is used
+
+// pml begin
+ppml[]+=pp[2];
+ppml[]+=newp; Point(newp) = {xb/cc+PmlDelta, 0, 0, lcair};
+ppml[]+=newp; Point(newp) = {xb/cc+PmlDelta, L, 0, lcair};
+ppml[]+=pp[3];
+
+For k In {0:2}
+ lpml1[]+=newl; Line(newl) = {ppml[k],ppml[k+1]};
+EndFor
+ppml2[]+=newp; Point(newp) = {0, xb/cc+L+PmlDelta, 0, lcair};
+
+lpml2[]+=newl; Circle(newl) = {ppml[2],p_[0],ppml2[0]};
+lpml2[]+=newl; Line(newl) = {ppml2[0],pp[4]};
+//pml end
+
+surfAirLL = newll ; Line Loop(surfAirLL) = {lbox[],-interfaceMonopoleAir[]};
+surfAir = news ; Plane Surface(surfAir) = {surfAirLL};
+
+// pml begin
+surfPmlLL1 = newll ; Line Loop(surfPmlLL1) = {lpml1[],-interfaceAirPML[0] };
+surfPml1 = news; Plane Surface(surfPml1) = {surfPmlLL1};
+
+surfPmlLL2 = newll ; Line Loop(surfPmlLL2) = {-lpml1[2],lpml2[],-interfaceAirPML[1] };
+surfPml2 = news; Plane Surface(surfPml2) = {surfPmlLL2};
+
+If(_use_structured_pml)
+ bndPml() = Boundary{Surface{surfPml};};
+
+ nn90 = Ceil[Pi*(rb+PmlDelta)/2/lcair];
+ nnDelta = Ceil[PmlDelta/lcair];
+
+ Transfinite Line {bndPml({0,2})} = nnDelta;
+ Transfinite Line {bndPml({1,3})} = nn90;
+ Transfinite Surface{surfPml};
+EndIf
+// pml_end
+
+
+// For aestetics
+Mesh.Light = 0; // so that to see mesh colors in the coax part (normal to surface has to be outward to see the colors with light :-))
+
+Recursive Color SkyBlue { Surface{surfAir};}
+Recursive Color Black { Surface{surfMonopole};}
+Recursive Color Yellow { Surface{surfPml1};}
+Recursive Color Yellow { Surface{surfPml2};}
+
+//=================================================
+// Physical regions
+//=================================================
+
+// surfaces
+Physical Surface(MONOPOLE) = {surfMonopole[]}; // MONOPOLE
+Physical Surface(AIR) = surfAir; // AIR
+Physical Surface(PML) = {surfPml1,surfPml2}; // PML
+
+// boundaries
+Physical Line(TERMINAL) = {feedMonople[]}; // Feeding
+Physical Line(CUT) = {lbox[{0}]};
+Physical Line(AXIS) = {lbox[{4}],axisMonopole[],lpml2[{1}]};
+
+// GROUND and SURFAIRINF
+If (Flag_AnalysisTypeExtBnd == 1) // SURFAIRINF is set to terminal ground (PEC type)
+ Physical Line(GROUND) = {lbox[{1}],lpml1[{0}],lpml1[{1}],lpml2[{0}]};
+Else
+ Physical Line(SURFAIRINF) = {lpml1[{1}],lpml2[{0}]};
+ Physical Line(GROUND) = {lbox[{1}],lpml1[{0}]};
+EndIf
+
diff --git a/ECE_full_wave_E-H_monopole/geo/monopole2D_cap_rad.geo b/ECE_full_wave_E-H_monopole/geo/monopole2D_cap_rad.geo
new file mode 100644
index 0000000000000000000000000000000000000000..8b83719405bc8ce32a70f2064216cc81712515d8
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/geo/monopole2D_cap_rad.geo
@@ -0,0 +1,107 @@
+Include "../monopole_data.pro";
+
+fac = s;
+flag = 1; // if 0, extruded mesh becomes free in the monopole
+
+lcdAntennaZero = a/4;
+lcdAntennaA = lcdAntennaZero;
+
+If (_use_fineMesh == 1)
+ lambdaMin = c0/fmax;
+ lcd = fac*lambdaMin/nopPerLambdaAir;
+ lcair = lcd;
+ nopL = flag*(Ceil[L/lcdAntennaZero]); // to make the antenna very fine along L
+Else
+ lcdAntennaZero = 4*lcdAntennaZero;
+ lambdaMin = c0/fmax;
+ lcair = fac*lambdaMin/nopPerLambdaAir;
+
+ delta = Sqrt(2.0*nuMonopole/(2*Pi*Freq*sigmaMonopole));
+ lcdAntennaA = fac*delta/nbDelta;
+
+ lcair1 = lambdaMin/10;
+ nopL = fac*6*flag*(Ceil[L/lcair1]);
+EndIf
+nopL = (nopL < 10) ? 10 : nopL;
+
+Printf("nopL %g",nopL);
+Printf("lcdAntennaA %g lcdAntennaZero %g ratio %g", lcdAntennaA, lcdAntennaZero, lcdAntennaZero/lcdAntennaA);
+
+//=================================================
+// Monopole
+//=================================================
+
+d = 0; // no gap
+yd = (d >= 0)? d:0.;
+p0 = newp ; Point(p0) = {0, yd, 0, lcdAntennaZero};
+p1 = newp ; Point(p1) = {a, yd, 0, lcdAntennaA};
+LineMonopoleBottom = newl ; Line(LineMonopoleBottom) = {p0,p1};
+
+surf[] = Extrude {0, L, 0} {
+ Line{LineMonopoleBottom}; Layers{nopL} ;
+};
+surfMonopole[0] = surf[1] ; LineMonopoleTop = surf[0] ;
+
+p_[] = Boundary{Line{LineMonopoleTop};};
+
+interfaceMonopoleAir[] = Boundary{ Surface{surfMonopole[]};} ;
+
+axisMonopole[] = interfaceMonopoleAir[{3}] ;
+feedMonople[] = interfaceMonopoleAir[{0}] ;
+interfaceMonopoleAir[] -= axisMonopole[];
+interfaceMonopoleAir[] -= feedMonople[];
+
+//=================================================
+// Air
+//=================================================
+
+pZero = newp; Point(pZero) = { 0, 0, 0, lcdAntennaZero};
+
+pp[]+=p1;
+pp[]+=newp; Point(newp) = { b, 0, 0, lcdAntennaZero};
+
+// reduce air around
+cc = 1;
+pp[]+=newp; Point(newp) = { xb/cc, 0, 0, lcair};
+pp[]+=newp; Point(newp) = { xb/cc, L, 0, lcair};
+pp[]+=newp; Point(newp) = { 0, xb/cc+L, 0, lcair};
+pp[]+=p_[0];
+
+For k In {0:2}
+ lbox[]+=newl; Line(newl) = {pp[k],pp[k+1]};
+EndFor
+lbox[]+=newl; Circle(newl) = {pp[3],p_[0],pp[4]};
+lbox[]+=newl; Line(newl) = {pp[4],pp[5]};
+interfaceAirPML[] = lbox[{2,3}]; // it will be the inf boundary if no PML is used
+
+surfAirLL = newll ; Line Loop(surfAirLL) = {lbox[],-interfaceMonopoleAir[]};
+surfAir = news ; Plane Surface(surfAir) = {surfAirLL};
+
+
+// For aestetics
+Mesh.Light = 0; // so that to see mesh colors in the coax part
+
+Recursive Color SkyBlue { Surface{surfAir};}
+Recursive Color Black { Surface{surfMonopole};}
+
+
+//=================================================
+// Physical regions
+//=================================================
+
+// surfaces
+Physical Surface(MONOPOLE) = {surfMonopole[]};
+Physical Surface(AIR) = surfAir;
+
+// boundaries
+Physical Line(TERMINAL) = {feedMonople[]};
+Physical Line(CUT) = {lbox[{0}]};
+
+Physical Line(AXIS) = {lbox[{4}],axisMonopole[]};
+
+If (Flag_AnalysisTypeExtBnd == 1)
+ Physical Line(GROUND) = {lbox[{1}],lbox[{2}],lbox[{3}]};
+Else
+ Physical Line(SURFAIRINF) = interfaceAirPML[];
+ Physical Line(GROUND) = {lbox[{1}]};
+EndIf
diff --git a/ECE_full_wave_E-H_monopole/geo/monopole2D_cyl_pml.geo b/ECE_full_wave_E-H_monopole/geo/monopole2D_cyl_pml.geo
new file mode 100644
index 0000000000000000000000000000000000000000..ed9874cb9d46f6ec8b53fdb87abb1e14afccaa53
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/geo/monopole2D_cyl_pml.geo
@@ -0,0 +1,144 @@
+Include "../monopole_data.pro";
+
+fac = s;
+flag = 1; // if 0, extruded mesh becomes free in the monopole
+_use_structured_pml = 0 ;
+
+lcdAntennaZero = a/4;
+lcdAntennaA = lcdAntennaZero;
+
+If (_use_fineMesh == 1)
+ lambdaMin = c0/fmax;
+ lcd = fac*lambdaMin/nopPerLambdaAir;
+ lcair = lcd;
+
+ nopL = flag*(Ceil[L/lcdAntennaZero]); // to make the antenna very fine along L
+Else
+ lcdAntennaZero = 4*lcdAntennaZero;
+ lambdaMin = c0/fmax;
+ lcair = fac*lambdaMin/nopPerLambdaAir;
+
+ delta = Sqrt(2.0*nuMonopole/(2*Pi*Freq*sigmaMonopole));
+ lcdAntennaA = fac*delta/nbDelta;
+
+ lcair1 = lambdaMin/10;
+ nopL = fac*6*flag*(Ceil[L/lcair1]);
+EndIf
+
+Printf("nopL %g",nopL);
+Printf("lcdAntennaA %g lcdAntennaZero %g ratio %g", lcdAntennaA, lcdAntennaZero, lcdAntennaZero/lcdAntennaA);
+
+//=================================================
+// Monopole
+//=================================================
+d = 0; // no gap
+yd = (d >= 0)? d:0.;
+p0 = newp ; Point(p0) = {0, yd, 0, lcdAntennaZero};
+p1 = newp ; Point(p1) = {a, yd, 0, lcdAntennaA};
+LineMonopoleBottom = newl ; Line(LineMonopoleBottom) = {p0,p1};
+
+surf[] = Extrude {0, L, 0} {
+ Line{LineMonopoleBottom}; Layers{nopL} ;
+};
+surfMonopole[0] = surf[1] ; LineMonopoleTop = surf[0] ;
+
+
+p_[] = Boundary{Line{LineMonopoleTop};};
+
+interfaceMonopoleAir[] = Boundary{ Surface{surfMonopole[]};} ;
+
+axisMonopole[] = interfaceMonopoleAir[{3}] ;
+feedMonople[] = interfaceMonopoleAir[{0}] ;
+interfaceMonopoleAir[] -= axisMonopole[];
+interfaceMonopoleAir[] -= feedMonople[];
+
+//=================================================
+// Air and PML
+//=================================================
+
+pZero = newp; Point(pZero) = { 0, 0, 0, lcdAntennaZero};
+
+pp[]+=p1;
+pp[]+=newp; Point(newp) = { b, 0, 0, lcdAntennaZero};
+
+pp[]+=newp; Point(newp) = { xb, 0, 0, lcair};
+pp[]+=newp; Point(newp) = { xb, yb, 0, lcair};
+pp[]+=newp; Point(newp) = { 0, yb, 0, lcair};
+pp[]+=p_[0];
+
+For k In {0:4}
+ lbox[]+=newl; Line(newl) = {pp[k],pp[k+1]};
+EndFor
+interfaceAirPML[] = lbox[{2,3}]; // it will be the inf boundary if no PML is used
+
+// pml begin
+ppml[]+=pp[2];
+ppml[]+=newp; Point(newp) = {xb+PmlDelta, 0, 0, lcair};
+ppml[]+=newp; Point(newp) = {xb+PmlDelta, yb+PmlDelta, 0, lcair};
+ppml[]+=pp[3];
+
+For k In {0:2}
+ lpml1[]+=newl; Line(newl) = {ppml[k],ppml[k+1]};
+EndFor
+
+ppml2[]+=newp; Point(newp) = {0, yb+PmlDelta, 0, lcair};
+lpml2[]+=newl; Line(newl) = {ppml[2],ppml2[0]};
+lpml2[]+=newl; Line(newl) = {ppml2[0],pp[4]};
+//pml end
+
+surfAirLL = newll ; Line Loop(surfAirLL) = {lbox[],-interfaceMonopoleAir[]};
+surfAir = news ; Plane Surface(surfAir) = {surfAirLL};
+
+
+// pml begin
+surfPmlLL1 = newll ; Line Loop(surfPmlLL1) = {lpml1[],-interfaceAirPML[0] };
+surfPml1 = news; Plane Surface(surfPml1) = {surfPmlLL1};
+
+surfPmlLL2 = newll ; Line Loop(surfPmlLL2) = {-lpml1[2],lpml2[],-interfaceAirPML[1] };
+surfPml2 = news; Plane Surface(surfPml2) = {surfPmlLL2};
+
+
+If(_use_structured_pml)
+ bndPml() = Boundary{Surface{surfPml};};
+
+ nn90 = Ceil[Pi*(rb+PmlDelta)/2/lcair];
+ nnDelta = Ceil[PmlDelta/lcair];
+
+ Transfinite Line {bndPml({0,2})} = nnDelta;
+ Transfinite Line {bndPml({1,3})} = nn90;
+ Transfinite Surface{surfPml};
+ // Recombine Surface {surfPml};
+EndIf
+// pml_end
+
+
+// For aestetics
+Mesh.Light = 0; // so that to see mesh colors in the coax part
+
+Recursive Color SkyBlue { Surface{surfAir};}
+Recursive Color Black { Surface{surfMonopole};}
+Recursive Color Yellow { Surface{surfPml1};}
+Recursive Color Yellow { Surface{surfPml2};}
+
+//=================================================
+// Physical regions
+//=================================================
+
+// surfaces
+Physical Surface(MONOPOLE) = {surfMonopole[]};
+Physical Surface(AIR) = surfAir;
+Physical Surface(PML) = {surfPml1,surfPml2};
+
+// boundaries
+Physical Line(TERMINAL) = {feedMonople[]}; // Feeding
+Physical Line(CUT) = {lbox[{0}]};
+
+Physical Line(AXIS) = {lbox[{4}],axisMonopole[],lpml2[{1}]};
+
+If (Flag_AnalysisTypeExtBnd == 1)
+ Physical Line(GROUND) = {lbox[{1}],lpml1[{0}],lpml1[{1}],lpml2[{0}]};
+Else
+ Physical Line(SURFAIRINF) = {lpml1[{1}],lpml2[{0}]};
+ Physical Line(GROUND) = {lbox[{1}],lpml1[{0}]};
+EndIf
+
diff --git a/ECE_full_wave_E-H_monopole/geo/monopole2D_cyl_rad.geo b/ECE_full_wave_E-H_monopole/geo/monopole2D_cyl_rad.geo
new file mode 100644
index 0000000000000000000000000000000000000000..f2f4b3ce149769e8a84d2124cec95c1a0b31bdb2
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/geo/monopole2D_cyl_rad.geo
@@ -0,0 +1,104 @@
+Include "../monopole_data.pro";
+
+fac = s;
+flag = 1; // if 0, extruded mesh becomes free in the monopole
+
+lcdAntennaZero = a/4;
+lcdAntennaA = lcdAntennaZero;
+
+If (_use_fineMesh == 1)
+ lambdaMin = c0/fmax;
+ lcd = fac*lambdaMin/nopPerLambdaAir;
+ lcair = lcd;
+
+ nopL = flag*(Ceil[L/lcdAntennaZero]); // to make the antenna very fine along L
+Else
+ lcdAntennaZero = 4*lcdAntennaZero;
+ lambdaMin = c0/fmax;
+ lcair = fac*lambdaMin/nopPerLambdaAir;
+
+ delta = Sqrt(2.0*nuMonopole/(2*Pi*Freq*sigmaMonopole));
+ lcdAntennaA = fac*delta/nbDelta;
+
+ lcair1 = lambdaMin/10;
+ nopL = fac*6*flag*(Ceil[L/lcair1]);
+EndIf
+nopL = (nopL < 10) ? 10 : nopL;
+
+Printf("nopL %g",nopL);
+Printf("lcdAntennaA %g lcdAntennaZero %g ratio %g", lcdAntennaA, lcdAntennaZero, lcdAntennaZero/lcdAntennaA);
+
+//=================================================
+// Monopole
+//=================================================
+
+d = 0; // no gap
+yd = (d >= 0)? d:0.;
+p0 = newp ; Point(p0) = {0, yd, 0, lcdAntennaZero};
+p1 = newp ; Point(p1) = {a, yd, 0, lcdAntennaA};
+LineMonopoleBottom = newl ; Line(LineMonopoleBottom) = {p0,p1};
+
+surf[] = Extrude {0, L, 0} {
+ Line{LineMonopoleBottom}; Layers{nopL} ;
+};
+surfMonopole[0] = surf[1] ; LineMonopoleTop = surf[0] ;
+
+p_[] = Boundary{Line{LineMonopoleTop};};
+
+interfaceMonopoleAir[] = Boundary{ Surface{surfMonopole[]};} ;
+
+axisMonopole[] = interfaceMonopoleAir[{3}] ;
+feedMonople[] = interfaceMonopoleAir[{0}] ;
+interfaceMonopoleAir[] -= axisMonopole[];
+interfaceMonopoleAir[] -= feedMonople[];
+
+//=================================================
+// Air
+//=================================================
+
+pZero = newp; Point(pZero) = { 0, 0, 0, lcdAntennaZero};
+
+pp[]+=p1;
+pp[]+=newp; Point(newp) = { b, 0, 0, lcdAntennaZero};
+
+pp[]+=newp; Point(newp) = { xb, 0, 0, lcair};
+pp[]+=newp; Point(newp) = { xb, yb, 0, lcair};
+pp[]+=newp; Point(newp) = { 0, yb, 0, lcair};
+pp[]+=p_[0];
+
+For k In {0:4}
+ lbox[]+=newl; Line(newl) = {pp[k],pp[k+1]};
+EndFor
+interfaceAirPML[] = lbox[{2,3}]; // it will be the inf boundary if no PML is used
+
+surfAirLL = newll ; Line Loop(surfAirLL) = {lbox[],-interfaceMonopoleAir[]};
+surfAir = news ; Plane Surface(surfAir) = {surfAirLL};
+
+
+// For aestetics
+Mesh.Light = 0; // so that to see mesh colors in the coax part
+
+Recursive Color SkyBlue { Surface{surfAir};}
+Recursive Color Black { Surface{surfMonopole};}
+
+
+//=================================================
+// Physical regions
+//=================================================
+
+// surfaces
+Physical Surface(MONOPOLE) = {surfMonopole[]};
+Physical Surface(AIR) = surfAir;
+
+// boundaries
+Physical Line(TERMINAL) = {feedMonople[]}; // Feeding
+Physical Line(CUT) = {lbox[{0}]};
+
+Physical Line(AXIS) = {lbox[{4}],axisMonopole[]};
+
+If (Flag_AnalysisTypeExtBnd == 1)
+ Physical Line(GROUND) = {lbox[{1}],lbox[{2}],lbox[{3}]};
+Else
+ Physical Line(SURFAIRINF) = interfaceAirPML[];
+ Physical Line(GROUND) = {lbox[{1}]};
+EndIf
diff --git a/ECE_full_wave_E-H_monopole/geo/monopole2D_sph_pml.geo b/ECE_full_wave_E-H_monopole/geo/monopole2D_sph_pml.geo
new file mode 100644
index 0000000000000000000000000000000000000000..dfcdda2554b09bcb15114ed130a98da5e3c53376
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/geo/monopole2D_sph_pml.geo
@@ -0,0 +1,124 @@
+Include "../monopole_data.pro";
+
+fac = s;
+flag = 1; // if 0, extruded mesh becomes free in the monopole
+_use_structured_pml = 1 ;
+
+If (_use_fineMesh == 1)
+ lcdAntennaZero = a/4;
+ lcdAntennaA = lcdAntennaZero;
+ lambdaMin = c0/fmax; // the same mesh for all the frequencies
+ lcd = fac*lambdaMin/nopPerLambdaAir;
+ lcair = lcd;
+
+ nopL = flag*(Ceil[L/lcdAntennaZero]);
+Else
+ lcdAntennaZero = a;
+ lambdaMin = lambda; // for a frequency dependent mesh, you should uncomment this line
+ lcair = fac*lambdaMin/nopPerLambdaAir;
+
+ no90 = Ceil[Pi*(rb)/2/lcair];
+ lcair = (no90 < 10)? Pi*(rb)/2/10 : no90;
+
+ delta = Sqrt(2.0*nuMonopole/(2*Pi*Freq*sigmaMonopole));
+ lcdAntennaA = fac*delta/nbDelta;
+
+ lcair1 = lambdaMin/(nopPerLambdaAir*2);
+ nopL = fac*3*flag*(Ceil[L/lcair1]);
+EndIf
+nopL = (nopL < 10) ? 10 : nopL;
+
+Printf("nopL %g",nopL);
+Printf("lcdAntennaA %g lcdAntennaZero %g ratio %g", lcdAntennaA, lcdAntennaZero, lcdAntennaZero/lcdAntennaA);
+
+
+//=================================================
+// Monopole
+//=================================================
+
+d = 0; // no gap
+yd = (d >= 0)? d:0.;
+p0 = newp ; Point(p0) = {0, yd, 0, lcdAntennaZero};
+p1 = newp ; Point(p1) = {a, yd, 0, lcdAntennaA};
+LineMonopoleBottom = newl ; Line(LineMonopoleBottom) = {p0,p1};
+
+surf[] = Extrude {0, L, 0} {
+ Line{LineMonopoleBottom}; Layers{nopL} ;
+};
+surfMonopole[0] = surf[1] ; LineMonopoleTop = surf[0] ;
+
+p_[] = Boundary{Line{LineMonopoleTop};};
+
+interfaceMonopoleAir[] = Boundary{ Surface{surfMonopole[]};} ;
+
+axisMonopole[] = interfaceMonopoleAir[{3}] ;
+feedMonople[] = interfaceMonopoleAir[{0}] ;
+interfaceMonopoleAir[] -= axisMonopole[];
+interfaceMonopoleAir[] -= feedMonople[];
+
+//=================================================
+// Air and PML
+//=================================================
+
+pZero = newp; Point(pZero) = { 0, 0, 0, lcdAntennaZero};
+
+pp[]+=p1;
+pp[]+=newp; Point(newp) = { b, 0, 0, lcdAntennaZero};
+
+pp[]+=newp; Point(newp) = { rb, 0, 0, lcair};
+pp[]+=newp; Point(newp) = { 0, rb, 0, lcair};
+pp[]+=p_[0];
+
+For k In {0:1}
+ lbox[]+=newl; Line(newl) = {pp[k],pp[k+1]};
+EndFor
+lbox[]+=newl; Circle(newl) = {pp[2],pZero,pp[3]};
+lbox[]+=newl; Line(newl) = {pp[3],pp[4]};
+interfaceAirPML[] = lbox[{2}]; // it will be the inf boundary if no PML is used
+
+// pml begin
+ppml[]+=pp[2];
+ppml[]+=newp; Point(newp) = {rb+PmlDelta, 0, 0, lcair};
+ppml[]+=newp; Point(newp) = {0, rb+PmlDelta, 0, lcair};
+ppml[]+=pp[3];
+
+lpml[]+=newl; Line(newl) = {ppml[0],ppml[1]};
+lpml[]+=newl; Circle(newl) = {ppml[1],pZero,ppml[2]};
+lpml[]+=newl; Line(newl) = {ppml[2],ppml[3]};
+//pml end
+
+surfAirLL = newll ; Line Loop(surfAirLL) = {lbox[],-interfaceMonopoleAir[]};
+surfAir = news ; Plane Surface(surfAir) = {surfAirLL};
+
+// pml begin
+surfPmlLL = newll ; Line Loop(surfPmlLL) = {lpml[],-interfaceAirPML[] };
+surfPml = news; Plane Surface(surfPml) = {surfPmlLL};
+
+If(_use_structured_pml)
+ bndPml() = Boundary{Surface{surfPml};};
+
+ nn90 = Ceil[Pi*(rb+PmlDelta)/2/lcair];
+ nnDelta = Ceil[PmlDelta/lcair];
+ Printf("*** *** *** ==================================================================> nnDelta = %e",nnDelta);
+ If (nnDelta < 4)
+ nnDelta = 4;
+ EndIf
+ Transfinite Line {bndPml({0,2})} = nnDelta;
+ Transfinite Line {bndPml({1,3})} = nn90;
+ Transfinite Surface{surfPml};
+ // Recombine Surface {surfPml}; // to obtain quadrangles - it does not work for formulation in E order 2 - because BF_Edge_4E is not ready for this
+EndIf
+// pml_end
+
+// For aestetics
+Mesh.Light = 0; // so that to see mesh colors in the coax part
+
+Recursive Color SkyBlue { Surface{surfAir};}
+Recursive Color Black { Surface{surfMonopole};}
+Recursive Color Yellow { Surface{surfPml};}
+
+//=================================================
+// Physical regions
+//=================================================
+
+Include "PhysicalRegions.geo";
diff --git a/ECE_full_wave_E-H_monopole/geo/monopole2D_sph_rad.geo b/ECE_full_wave_E-H_monopole/geo/monopole2D_sph_rad.geo
new file mode 100644
index 0000000000000000000000000000000000000000..31ac3a1abc5b8516d2477ce1882cd66c14373f8d
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/geo/monopole2D_sph_rad.geo
@@ -0,0 +1,82 @@
+Include "../monopole_data.pro";
+
+fac = s;
+flag = 1; // if 0, extruded mesh becomes free in the monopole
+
+If (_use_fineMesh == 1)
+ lcdAntennaZero = a/4;
+ lcdAntennaA = lcdAntennaZero;
+ lambdaMin = c0/fmax; // the same mesh for all the frequencies
+ lcd = fac*lambdaMin/nopPerLambdaAir;
+ lcair = lcd;
+
+ nopL = flag*(Ceil[L/lcdAntennaZero]); // to make the antenna very fine along L
+Else
+ lcdAntennaZero = a;
+ lambdaMin = lambda; // for a frequency dependent mesh, you should uncomment this line
+ lcair = fac*lambdaMin/nopPerLambdaAir;
+
+ no90 = Ceil[Pi*(rb)/2/lcair];
+ lcair = (no90 < 10)? Pi*(rb)/2/10 : no90;
+
+ delta = Sqrt(2.0*nuMonopole/(2*Pi*Freq*sigmaMonopole));
+ lcdAntennaA = fac*delta/nbDelta;
+
+ lcair1 = lambdaMin/(nopPerLambdaAir*2);
+ nopL = fac*3*flag*(Ceil[L/lcair1]);
+EndIf
+nopL = (nopL < 10) ? 10 : nopL;
+
+Printf("====================================> nopL %g",nopL);
+Printf("lcdAntennaA %g lcdAntennaZero %g ratio %g", lcdAntennaA, lcdAntennaZero, lcdAntennaZero/lcdAntennaA);
+
+//=================================================
+// Monopole - no gap
+//=================================================
+
+p0 = newp ; Point(p0) = {0, 0, 0, lcdAntennaZero};
+p1 = newp ; Point(p1) = {a, 0, 0, lcdAntennaA};
+LineMonopoleBottom = newl ; Line(LineMonopoleBottom) = {p0,p1};
+
+surf[] = Extrude {0, L, 0} {
+ Line{LineMonopoleBottom}; Layers{nopL} ;
+};
+surfMonopole[0] = surf[1] ; LineMonopoleTop = surf[0] ;
+
+p_[] = Boundary{Line{LineMonopoleTop};};
+
+interfaceMonopoleAir[] = Boundary{ Surface{surfMonopole[]};} ;
+
+axisMonopole[] = interfaceMonopoleAir[{3}] ;
+feedMonople[] = interfaceMonopoleAir[{0}] ;
+interfaceMonopoleAir[] -= axisMonopole[];
+interfaceMonopoleAir[] -= feedMonople[];
+
+//=================================================
+// Air
+//=================================================
+
+pZero = newp; Point(pZero) = { 0, 0, 0, lcdAntennaZero};
+
+pp[]+=p1;
+pp[]+=newp; Point(newp) = { b, 0, 0, lcdAntennaZero};
+
+pp[]+=newp; Point(newp) = { rb, 0, 0, lcair};
+pp[]+=newp; Point(newp) = { 0, rb, 0, lcair};
+pp[]+=p_[0];
+
+For k In {0:1}
+ lbox[]+=newl; Line(newl) = {pp[k],pp[k+1]};
+EndFor
+lbox[]+=newl; Circle(newl) = {pp[2],pZero,pp[3]};
+lbox[]+=newl; Line(newl) = {pp[3],pp[4]};
+interfaceAirPML[] = lbox[{2}]; // it will be the inf boundary if no PML is used
+
+surfAirLL = newll ; Line Loop(surfAirLL) = {lbox[],-interfaceMonopoleAir[]};
+surfAir = news ; Plane Surface(surfAir) = {surfAirLL};
+
+
+//=================================================
+// Physical regions
+//=================================================
+Include "PhysicalRegions.geo";
diff --git a/ECE_full_wave_E-H_monopole/geo/monopole2Dcoax_sph_rad.geo b/ECE_full_wave_E-H_monopole/geo/monopole2Dcoax_sph_rad.geo
new file mode 100644
index 0000000000000000000000000000000000000000..c4417fe64c7b987fddbd142e1e7d34f0d350e63b
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/geo/monopole2Dcoax_sph_rad.geo
@@ -0,0 +1,131 @@
+Include "../monopole_data.pro";
+
+fac = s;
+
+// characteristic lengths & some transfinite number of divisions
+nopPerLambda = NOP_perLAMBDA_AIR;
+lc = lambda/nopPerLambda ; // rule of thumb (min 10 divisions)
+
+lcd = (lc < a/2) ? lc :a/2 ;
+flag = 1; // if 0, extruded mesh becomes free
+
+lcdAntennaZero = a/4;
+lcdAntennaA = lcdAntennaZero;
+
+nopL = flag*Ceil[L/lcd/2/fac] ;
+nopdCoax = flag*Ceil[Abs[d]*nopL/L] ;
+
+If (_use_fineMesh == 1)
+ lambdaMin = c0/fmax;
+ lcd = fac*lambdaMin/nopPerLambdaAir;
+ lcair = lcd;
+
+ nopL = flag*(Ceil[L/lcdAntennaZero]); // to make the antenna very fine along L
+ nopdCoax = flag*Ceil[Abs[d]*nopL/L] ;
+Else
+ lambdaMin = c0/fmax;
+ lcair = fac*lambdaMin/nopPerLambdaAir;
+ lcdAntennaZero = lcair;
+ xxx = 1;
+ lcdAntennaZero = lcdAntennaZero/xxx;
+ lcdAntennaA = lcair;
+
+ lcair1 = lcair;
+ lcair1 = lcair1/xxx;
+ nopL = fac*flag*(Ceil[L/lcair1]);
+ nopdCoax = flag*Ceil[Abs[d]*nopL/L] ;
+EndIf
+
+Printf("nopL %g",nopL);
+Printf("lcdAntennaA %g lcdAntennaZero %g ratio %g", lcdAntennaA, lcdAntennaZero,lcdAntennaA/lcdAntennaZero);
+
+lcdAntennaZeroB = lcdAntennaZero;
+
+//=================================================
+// Monopole
+//=================================================
+
+p0 = newp ; Point(p0) = {0, 0, 0, lcdAntennaZero};
+p1 = newp ; Point(p1) = {a, 0, 0, lcdAntennaA};
+
+LineMonopoleBottom = newl ; Line(LineMonopoleBottom) = {p0,p1};
+
+surf[] = Extrude {0, L, 0} {
+ Line{LineMonopoleBottom}; Layers{nopL} ;
+};
+surfMonopole[0] = surf[1] ; LineMonopoleTop = surf[0] ;
+
+p_[] = Boundary{Line{LineMonopoleTop};};
+
+interfaceMonopoleAir[] = Boundary{ Surface{surfMonopole[]};} ;
+
+axisMonopole[] = interfaceMonopoleAir[{3}] ;
+feedMonople[] = interfaceMonopoleAir[{0}] ;
+interfaceMonopoleAir[] -= axisMonopole[];
+interfaceMonopoleAir[] -= feedMonople[];
+
+//=================================================
+// Air, circular bnd, no pml
+//=================================================
+
+pZero = newp; Point(pZero) = { 0, 0, 0, lcd};
+
+pp[]+=p1;
+pp[]+=newp; Point(newp) = { b, 0, 0, lcdAntennaA};
+pp[]+=newp; Point(newp) = { rb, 0, 0, lcair};
+pp[]+=newp; Point(newp) = { 0, rb, 0, lcair};
+pp[]+=p_[0];
+
+For k In {0:1}
+ lbox[]+=newl; Line(newl) = {pp[k],pp[k+1]};
+EndFor
+lbox[]+=newl; Circle(newl) = {pp[2],pZero,pp[3]};
+lbox[]+=newl; Line(newl) = {pp[3],pp[4]};
+interfaceAirPML[] = lbox[{2}];
+
+//=================================================
+// add the COAX part
+//=================================================
+// the conductor
+surfCoax[] = Extrude {0, d, 0} {
+ Line{LineMonopoleBottom}; Layers{nopdCoax} ;
+};
+surfCoaxConductor[0] = surfCoax[1] ; feedMonople[] = surfCoax[0] ;
+surfCoaxDiel[] = Extrude {0, d, 0} {
+ Line{lbox[0]}; Layers{nopdCoax} ;
+ };
+surfCoaxDielectric[0] = surfCoaxDiel[1] ; cut[] = surfCoaxDiel[0] ;
+
+surfAirLL = newll ; Line Loop(surfAirLL) = {lbox[],-interfaceMonopoleAir[]};
+surfAir = news ; Plane Surface(surfAir) = {surfAirLL};
+
+// For aestetics
+Mesh.Light = 0; // so that to see mesh colors in the coax part
+
+Recursive Color SkyBlue { Surface{surfAir};}
+Recursive Color Black { Surface{surfMonopole};}
+Recursive Color Orange { Surface{surfCoaxConductor[]};}
+Recursive Color Blue { Surface{surfCoaxDielectric[]};}
+
+//=================================================
+// Physical regions
+//=================================================
+
+// surfaces
+Physical Surface(MONOPOLE) = {surfMonopole[]};
+Physical Surface(AIR) = surfAir;
+Physical Surface(COAXCONDUCTOR) = {surfCoaxConductor[]};
+Physical Surface(COAXDIELECTRIC) = {surfCoaxDielectric};
+
+// boundaries
+Physical Line(TERMINAL) = {feedMonople[]}; // Feeding
+
+Physical Line(CUT) = {surfCoaxDiel[0]};
+Physical Line(AXIS) = {lbox[{3}],axisMonopole[],surfCoax[3]};
+
+If (Flag_AnalysisTypeExtBnd==1)
+ Physical Line(GROUND) = {surfCoaxDiel[2],lbox[{1}],lbox[{2}]};
+Else
+ Physical Line(SURFAIRINF) = interfaceAirPML;
+ Physical Line(GROUND) = {surfCoaxDiel[2],lbox[{1}]};
+EndIf
\ No newline at end of file
diff --git a/ECE_full_wave_E-H_monopole/geo/monopole2Dgap_sph_rad.geo b/ECE_full_wave_E-H_monopole/geo/monopole2Dgap_sph_rad.geo
new file mode 100644
index 0000000000000000000000000000000000000000..4d9e454fab57075047af6c2a41cee4d4eedb257d
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/geo/monopole2Dgap_sph_rad.geo
@@ -0,0 +1,113 @@
+Include "../monopole_data.pro";
+
+fac = s;
+nopPerLambda = NOP_perLAMBDA_AIR;
+lc = lambda/nopPerLambda ; // rule of thumb (min 10 divisions)
+lcd = (lc < a/2) ? lc :a/2 ;
+
+flag = 1; // if 0, extruded mesh becomes free
+lcdAntennaZero = a/4;
+lcdAntennaA = lcdAntennaZero;
+
+nopL = flag*Ceil[L/lcd/2/fac] ;
+
+If (_use_fineMesh == 1)
+ lambdaMin = c0/fmax;
+ lcd = fac*lambdaMin/nopPerLambdaAir;
+ lcair = lcd;
+
+ nopL = flag*(Ceil[L/lcdAntennaZero]); // to make the antenna very fine along L
+Else
+
+ lambdaMin = c0/fmax;
+ lcair = fac*lambdaMin/nopPerLambdaAir;
+ lcdAntennaZero = lcair;
+ xxx = 1;
+ lcdAntennaZero = lcdAntennaZero/xxx;
+ lcdAntennaA = lcair;
+
+ lcair1 = lcair;
+ lcair1 = lcair1/xxx;
+ nopL = fac*flag*(Ceil[L/lcair1]);
+ If (d < 0)
+ nopdCoax = flag*Ceil[Abs[d]*nopL/L] ;
+ EndIf
+EndIf
+
+Printf("====================================> nopL %g",nopL);
+Printf("lcdAntennaA %g lcdAntennaZero %g ratio %g", lcdAntennaA, lcdAntennaZero, lcdAntennaZero/lcdAntennaA);
+
+lcdAntennaZeroB = lcdAntennaZero;
+
+//=================================================
+// Monopole - gap > 0
+//=================================================
+
+// d is assumed > 0
+p0 = newp ; Point(p0) = {0, d, 0, lcdAntennaZero};
+p1 = newp ; Point(p1) = {a, d, 0, lcdAntennaZeroB};
+
+LineMonopoleBottom = newl ; Line(LineMonopoleBottom) = {p0,p1};
+
+surf[] = Extrude {0, L, 0} {
+ Line{LineMonopoleBottom}; Layers{nopL} ;
+};
+surfMonopole[0] = surf[1] ; LineMonopoleTop = surf[0] ;
+
+p_[] = Boundary{Line{LineMonopoleTop};};
+
+interfaceMonopoleAir[] = Boundary{ Surface{surfMonopole[]};} ;
+
+axisMonopole[] = interfaceMonopoleAir[{3}] ;
+feedMonople[] = interfaceMonopoleAir[{0}] ;
+interfaceMonopoleAir[] -= axisMonopole[];
+interfaceMonopoleAir[] -= feedMonople[];
+
+//=================================================
+// Air, circular truncation boundary
+//=================================================
+
+pZero = newp; Point(pZero) = { 0, 0, 0, lcd};
+
+pp[]+=p1;
+pp[]+=newp; Point(newp) = { a, 0, 0, lcdAntennaA};
+pp[]+=newp; Point(newp) = { rb, 0, 0, lcair};
+pp[]+=newp; Point(newp) = { 0, rb, 0, lcair};
+pp[]+=p_[0];
+
+For k In {0:1}
+ lbox[]+=newl; Line(newl) = {pp[k],pp[k+1]};
+EndFor
+lbox[]+=newl; Circle(newl) = {pp[2],pZero,pp[3]};
+lbox[]+=newl; Line(newl) = {pp[3],pp[4]};
+interfaceAirPML[] = lbox[{2}]; // it will be the inf boundary if no PML is used
+
+surfAirLL = newll ; Line Loop(surfAirLL) = {lbox[],-interfaceMonopoleAir[]};
+surfAir = news ; Plane Surface(surfAir) = {surfAirLL};
+
+// For aestetics
+Mesh.Light = 0; // so that to see mesh colors in the coax part
+
+Recursive Color SkyBlue { Surface{surfAir};}
+Recursive Color Black { Surface{surfMonopole};}
+
+//=================================================
+// Physical regions
+//=================================================
+
+// surfaces
+Physical Surface(MONOPOLE) = {surfMonopole[]};
+Physical Surface(AIR) = surfAir;
+
+// boundaries
+Physical Line(TERMINAL) = {feedMonople[]}; // Feeding
+Physical Line(CUT) = {lbox[{0}]};
+
+Physical Line(AXIS) = {lbox[{3}],axisMonopole[]};
+
+If (Flag_AnalysisTypeExtBnd==1)
+ Physical Line(GROUND) = {lbox[{1}],lbox[{2}]};
+Else
+ Physical Line(SURFAIRINF) = interfaceAirPML;
+ Physical Line(GROUND) = {lbox[{1}]};
+EndIf
diff --git a/ECE_full_wave_E-H_monopole/geo/monopole3D_sph_pml.geo b/ECE_full_wave_E-H_monopole/geo/monopole3D_sph_pml.geo
new file mode 100644
index 0000000000000000000000000000000000000000..14a0e1bd6a0a92c4da4801823f23f5cdd903d48a
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/geo/monopole3D_sph_pml.geo
@@ -0,0 +1,132 @@
+Include "../monopole_data.pro";
+
+//****************************************************** start of geometry
+SetFactory("OpenCASCADE");
+
+//Mesh.MeshSizeFactor = s ;
+// Mesh.Algorithm3D = 10; // 3D mesh algorithm (1=Delaunay, 4=Frontal, 5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D, 9=R-tree, 10=HXT)
+Mesh.Optimize = 1 ;
+
+ang_section = ang_section_degrees*Pi/180; // => cuts not be added if a section is considered
+
+idxCyl = newv; Cylinder(newv) = {0,0,0,0,0,L,a,ang_section}; // monopole along Oz
+idxDisk = news; Disk(news) = {0,0,0,b};
+idxSph = newv; Sphere(newv) = {0,0,0,rb,0,Pi/2,ang_section};
+idxSphPml = newv; Sphere(newv) = {0,0,0,rb+PmlDelta,0,Pi/2,ang_section};
+
+// For keeping the same rotation axis as in 2D
+Rotate{ {1, 0, 0}, {0, 0, 0}, -Pi/2 } { Surface{:};Volume{:}; }
+
+BooleanFragments{Volume{:}; Surface{:}; Delete;}{} // to intersect everything
+// Avoid the use of BooleanDifference, it is not trustworthy
+If(ang_section_degrees<360)
+ Recursive Delete{Surface{17};}
+EndIf
+
+vol() = Volume{:}; // Not really necessary but...
+Printf("Volumes in the geometry after BooleanFragments for volumes",vol());
+
+// Hunting the volume numbers from the interface (GUI: Tools -> Visibility -> list, click on volumes and see which is which
+volMonopole[] = vol(0);
+volAir[] = vol(1);
+volPml[] = vol(2);
+
+surf_bnd() = Abs(CombinedBoundary{Volume{:};});
+Printf("all boundary surfaces", surf_bnd());
+
+// to better see the surfaces,Tools -> Options -> Visibility -> Geometry -> Surfaces, Surface Labels and you will see the number of surfaces
+// to better see the selected surfaces, Tools -> Options -> Geometry -> Aspect -> Surface Display -> Solid
+
+If(ang_section_degrees == 360)
+ surfaceFeed[] = surf_bnd(0);
+ surfaceAirInf[] = surf_bnd(3);
+ surfaceGround[] = surf_bnd({1,4});
+ surfaceNonTerm[] = surf_bnd(2);
+EndIf
+If(ang_section_degrees < 360)
+ surfaceFeed[] = surf_bnd(0);
+ surfaceAirInf[] = surf_bnd(7);
+ surfaceGround[] = surf_bnd({4,9});
+ surfaceNonTerm[] = surf_bnd(6);
+ surfaceCutA[] = surf_bnd({1,5,10});
+ surfaceCutB[] = surf_bnd({2,3,8});
+EndIf
+
+//=================================================
+// Physical regions
+//=================================================
+// volumes
+Physical Volume("Monopole", MONOPOLE) = {volMonopole[]};
+Physical Volume("Air",AIR) = {volAir[]};
+Physical Volume("PML", PML) = {volPml[]};
+
+// surfaces
+Physical Surface("Terminal", TERMINAL) = {surfaceFeed[]};
+Physical Surface("Fictitious bnd - radiant", SURFAIRINF) = {surfaceAirInf[]};
+Physical Surface("Ground",GROUND) = {surfaceGround[]};
+Physical Surface("Non terminal", NONTERM) = {surfaceNonTerm[]};
+
+If(ang_section_degrees < 360)
+ // Most likely not needed for the E-formulation
+ Periodic Surface {surfaceCutA[]} = {surfaceCutB[]} Rotate {{0,1,0},{0,0,0}, -ang_section};
+
+ Physical Surface("Cut A", SURFAIRINF*10) = {surfaceCutA[]};
+ Physical Surface("Cut B", SURFAIRINF*10+1) = {surfaceCutB[]};
+EndIf
+
+If (_use_fineMesh == 1)
+ // Adapting the mesh if necesary
+ MeshSize{PointsOf{Volume{:};}} = L/10;
+ MeshSize{PointsOf{Volume{volMonopole[]};}} = L/40;
+ MeshSize{PointsOf{Volume{volPml[]};}} = s*lambda/nopPerLambdaAir/4;
+ MeshSize{PointsOf{Surface{surfaceAirInf[]};}} = s*lambda/nopPerLambdaAir/2;
+
+ // We can activate the calculation of mesh element sizes based on curvature:
+ // Mesh.MeshSizeFromCurvature = 1;
+
+ // And we set the minimum number of elements per 2*Pi radians:
+ Mesh.MinimumElementsPerTwoPi = (ang_section_degrees < 360) ? 4*360/ang_section_degrees : 24/2;
+ // Mesh.MeshSizeMax = s*lambda/nopPerLambdaAir/5;
+ // Mesh.MeshSizeMin = Mesh.MeshSizeMax/10;
+Else
+ lcair = s*lambda/nopPerLambdaAir;
+ nopL = (Ceil[L/(a/4)]); // to make the antenna very fine along L
+
+ no90 = Ceil[Pi*(rb)/2/lcair];
+ If (no90 < 10)
+ lcair = Pi*(rb)/2/10;
+ EndIf
+
+ MeshSize{PointsOf{Volume{:};}} = lcair; // does this mean that at first all the volumes have this setting? Why did you chose L/10 and not something depeding on the frequency?
+ Printf("L = %g",L);
+ Printf("L/10 = %g",L/10);
+ Printf("lcair = %g",lcair);
+ Printf("nopL = %g",nopL);
+ lcMonopol = L/10; // changed to something smaller so that I can solve it quicker on my laptop
+ MeshSize{PointsOf{Volume{volMonopole[]};}} = lcMonopol;
+ MeshSize{PointsOf{Volume{volPml[]};}} = lcair;
+ MeshSize{PointsOf{Surface{surfaceAirInf[]};}} = lcair;
+EndIf
+
+// Transfinite PML
+If(ang_section_degrees < 360)
+ gg = ang_section_degrees/30;
+ ff = 1;
+ Transfinite Line {21,22,23} = 4*ff; // width
+ Transfinite Line {2, 4, 18, 20} = 10*ff; // arcs radial
+ Transfinite Line {3,19} = 5*gg*ff;
+
+ Transfinite Surface {14,15,16}; // sides
+ Transfinite Surface {6} = {1,2,3};
+ Transfinite Surface {13} ={11,12,13};
+ Transfinite Volume {volPml[]};
+ _use_recombine=0;
+ If(_use_recombine)
+ Recombine Surface {14,15,16};
+ Recombine Surface {6,13};
+ Recombine Volume {volPml[]};
+ EndIf
+EndIf
+
+Cohomology(1) {{NONTERM},{}}; // this will generate the cut needed by the formulation in H, with no NONTERM + 1
+
diff --git a/ECE_full_wave_E-H_monopole/geo/monopole3D_sph_rad.geo b/ECE_full_wave_E-H_monopole/geo/monopole3D_sph_rad.geo
new file mode 100644
index 0000000000000000000000000000000000000000..2cb84ee0b811318967f0e6d2dad119d4e5ca316e
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/geo/monopole3D_sph_rad.geo
@@ -0,0 +1,115 @@
+Include "../monopole_data.pro";
+
+//****************************************************** start of geometry
+SetFactory("OpenCASCADE");
+
+//Mesh.MeshSizeFactor = s ;
+// Mesh.Algorithm3D = 10; // 3D mesh algorithm (1=Delaunay, 4=Frontal, 5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D, 9=R-tree, 10=HXT)
+Mesh.Optimize = 1 ;
+
+ang_section = ang_section_degrees*Pi/180; // => cuts not be added if a section is considered
+
+idxCyl = newv; Cylinder(idxCyl) = {0,0,0,0,0,L,a,ang_section}; // monopole along Oz
+idxDisk = news; Disk(idxDisk) = {0,0,0,b};
+idxSph = newv; Sphere(idxSph) = {0,0,0,rb,0,Pi/2,ang_section};
+
+// For keeping the same rotation axis as in 2D
+Rotate{ {1, 0, 0}, {0, 0, 0}, -Pi/2 } { Surface{:};Volume{:}; }
+
+BooleanFragments{Volume{:}; Surface{:}; Delete;}{} // to intersect everything
+// Avoid the use of BooleanDifference, it is not trustworthy
+If(ang_section_degrees<360)
+ Recursive Delete{Surface{13};}
+EndIf
+
+vol() = Volume{:}; // Not really necessary but...
+Printf("Volumes in the geometry after BooleanFragments for volumes",vol());
+
+// Hunting the volume numbers from the interface (GUI: Tools -> Visibility -> list, click on volumes and see which is which
+volMonopole[] = vol(0);
+volAir[] = vol(1);
+
+surf_bnd() = Abs(CombinedBoundary{Volume{:};});
+Printf("all boundary surfaces", surf_bnd());
+
+// to better see the surfaces,Tools -> Options -> Visibility -> Geometry -> Surfaces, Surface Labels and you will see the number of surfaces
+// to better see the selected surfaces, Tools -> Options -> Geometry -> Aspect -> Surface Display -> Solid
+
+If(ang_section_degrees == 360)
+ surfaceFeed[] = surf_bnd(0);
+ surfaceAirInf[] = surf_bnd(1);
+ surfaceGround[] = surf_bnd(2);
+ surfaceNonTerm[] = surf_bnd(3);
+EndIf
+If(ang_section_degrees < 360)
+ surfaceFeed[] = surf_bnd(0);
+ surfaceAirInf[] = surf_bnd(3);
+ surfaceGround[] = surf_bnd(5);
+ surfaceNonTerm[] = surf_bnd(7);
+ surfaceCutA[] = surf_bnd({1,6});
+ surfaceCutB[] = surf_bnd({2,4});
+EndIf
+
+//=================================================
+// Physical regions
+//=================================================
+// volumes
+Physical Volume("Monopole", MONOPOLE) = {volMonopole[]}; // MONOPOLE
+Physical Volume("Air",AIR) = {volAir[]}; // AIR
+// surfaces
+Physical Surface("Terminal", TERMINAL) = {surfaceFeed[]}; // FEED
+Physical Surface("Fictitious bnd - radiant", SURFAIRINF) = {surfaceAirInf[]};
+Physical Surface("Ground",GROUND) = {surfaceGround[]}; // GROUND
+Physical Surface("Non terminal", NONTERM) = {surfaceNonTerm[]}; // NON Terminal
+
+If(ang_section_degrees < 360)
+ // Most likely not needed for the E-formulation
+ Periodic Surface {surfaceCutA[]} = {surfaceCutB[]} Rotate {{0,1,0},{0,0,0}, -ang_section};
+ Physical Surface("Cut A", SURFAIRINF*10) = {surfaceCutA[]};
+ Physical Surface("Cut B", SURFAIRINF*10+1) = {surfaceCutB[]};
+EndIf
+
+If (_use_fineMesh == 1)
+ // Adapting the mesh if necesary
+ MeshSize{PointsOf{Volume{:};}} = L/10;
+ MeshSize{PointsOf{Volume{volMonopole[]};}} = L/40;
+ MeshSize{PointsOf{Surface{surfaceAirInf[]};}} = s*lambda/nopPerLambdaAir/2;
+ // We can activate the calculation of mesh element sizes based on curvature:
+ // Mesh.MeshSizeFromCurvature = 1;
+
+ // And we set the minimum number of elements per 2*Pi radians:
+ Mesh.MinimumElementsPerTwoPi = (ang_section_degrees < 360) ? 4*360/ang_section_degrees : 24/2;
+ // Mesh.MeshSizeMax = s*lambda/nopPerLambdaAir/5;
+ // Mesh.MeshSizeMin = Mesh.MeshSizeMax/10;
+Else
+ If(1) // frequency dependent lcair
+ lcair = s*lambda/nopPerLambdaAir;
+ Else
+ // the same mesh for all the freq
+ lambdaMin = c0/fmax;
+ lcair = s*lambdaMin/nopPerLambdaAir; // the same mesh for all the frequencies, not reccomended in 3d , especially for order 2
+ EndIf
+
+ //
+ nopL = (Ceil[L/(a/4)]); // to make the antenna very fine along L
+
+ no90 = Ceil[Pi*(rb)/2/lcair];
+ If (no90 < 10)
+ lcair = Pi*(rb)/2/10;
+ EndIf
+
+ MeshSize{PointsOf{Volume{:};}} = lcair; // does this mean that at first all the volumes have this setting? Why did you chose L/10 and not something depeding on the frequency?
+ Printf("L = %g",L);
+ Printf("L/10 = %g",L/10);
+ Printf("lcair = %g",lcair);
+ Printf("nopL = %g",nopL);
+ lcMonopol = L/10; // changed to something smaller so that I can solve it quicker on my laptop
+ MeshSize{PointsOf{Volume{volMonopole[]};}} = lcMonopol;
+ MeshSize{PointsOf{Surface{surfaceAirInf[]};}} = lcair;
+EndIf
+
+
+
+// cut for the H formulation
+
+Cohomology(1) {{NONTERM},{}}; // this will generate the cut needed by the formulation in H, with no NONTERM + 1
diff --git a/ECE_full_wave_E-H_monopole/geo/monopole3D_sph_rad_new.geo b/ECE_full_wave_E-H_monopole/geo/monopole3D_sph_rad_new.geo
new file mode 100644
index 0000000000000000000000000000000000000000..04193d8159e0d3e0c85a4ebd8f2cc5c6d6570b66
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/geo/monopole3D_sph_rad_new.geo
@@ -0,0 +1,103 @@
+Include "../monopole_data.pro";
+
+//****************************************************** start of geometry
+SetFactory("OpenCASCADE");
+
+//Mesh.MeshSizeFactor = s ;
+// Mesh.Algorithm3D = 10; // 3D mesh algorithm (1=Delaunay, 4=Frontal, 5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D, 9=R-tree, 10=HXT)
+Mesh.Optimize = 1 ;
+
+// And we set the minimum number of elements per 2*Pi radians:
+// We can activate the calculation of mesh element sizes based on curvature:
+Mesh.MeshSizeFromCurvature = 1;
+Mesh.MinimumElementsPerTwoPi = (_use_fineMesh ? 16 : 8) ; // 24
+
+// Mesh.MeshSizeMax = s*lambda/nopPerLambdaAir/5;
+// Mesh.MeshSizeMin = Mesh.MeshSizeMax/10;
+
+ang_section = ang_section_degrees*Pi/180; // => cuts not be added if a section is considered
+
+idxCyl = newv; Cylinder(idxCyl) = {0,0,0,0,0,L,a,ang_section}; // monopole along Oz
+// idxDisk = news; Disk(idxDisk) = {0,0,0,b};
+idxCylDisk = newv; Cylinder(idxCylDisk) = {0,0,0,0,0,L,b,ang_section}; // for better controlling the mesh
+idxSph = newv; Sphere(idxSph) = {0,0,0,rb,0,Pi/2,ang_section};
+
+// For keeping the same rotation axis as in 2D
+Rotate{ {1, 0, 0}, {0, 0, 0}, -Pi/2 } { Surface{:};Volume{:}; }
+
+BooleanFragments{Volume{:}; Surface{:}; Delete;}{} // to intersect everything
+
+vol() = Volume{:}; // Not really necessary but...
+Printf("Volumes in the geometry after BooleanFragments for volumes",vol());
+
+// Hunting the volume numbers from the interface (GUI: Tools -> Visibility -> list, click on volumes and see which is which
+volMonopole() = vol(0);
+volAir() = vol({1,2});
+
+surf_bnd() = Abs(CombinedBoundary{Volume{:};});
+Printf("all boundary surfaces", surf_bnd());
+
+If(ang_section_degrees == 360)
+ surfaceFeed() = surf_bnd(0);
+ surfaceAirInf() = surf_bnd(2);
+ surfaceGround() = surf_bnd(3);
+ surfaceNonTerm() = surf_bnd(1);
+EndIf
+If(ang_section_degrees < 360)
+ surfaceFeed() = surf_bnd(0);
+ surfaceAirInf() = surf_bnd(6);
+ surfaceGround() = surf_bnd(8);
+ surfaceNonTerm() = surf_bnd(4);
+ surfaceCutA() = surf_bnd({1,5,9});
+ surfaceCutB() = surf_bnd({2,3,7});
+EndIf
+
+//=================================================
+// Physical regions
+//=================================================
+// volumes
+Physical Volume("Monopole", MONOPOLE) = {volMonopole[]}; // MONOPOLE
+Physical Volume("Air",AIR) = {volAir[]}; // AIR
+// surfaces
+Physical Surface("Terminal", TERMINAL) = {surfaceFeed[]}; // FEED
+Physical Surface("Fictitious bnd - radiant", SURFAIRINF) = {surfaceAirInf[]};
+Physical Surface("Ground",GROUND) = {surfaceGround[]}; // GROUND
+Physical Surface("Non terminal", NONTERM) = {surfaceNonTerm[]}; // NON Terminal
+
+If(ang_section_degrees < 360)
+ // Most likely not needed for the E-formulation
+ Periodic Surface {surfaceCutA[]} = {surfaceCutB[]} Rotate {{0,1,0},{0,0,0}, -ang_section};
+ Physical Surface("Cut A", SURFAIRINF*10) = {surfaceCutA[]};
+ Physical Surface("Cut B", SURFAIRINF*10+1) = {surfaceCutB[]};
+EndIf
+
+
+// Adapting the mesh...
+// ==============================
+If(1) // frequency dependent lcair
+ lcair = s*lambda/nopPerLambdaAir;
+Else
+ // the same mesh for all the freq
+ lambdaMin = c0/fmax;
+ lcair = s*lambdaMin/nopPerLambdaAir; // the same mesh for all the frequencies, not reccomended in 3d , especially for order 2
+EndIf
+
+// nopL = Ceil[L/(a/4)]; // to make the antenna very fine along L
+lcMonopol = s*L/10;
+MeshSize{PointsOf{Volume{:};}} = (_use_fineMesh) ? lcMonopol*6 : lcMonopol*8;
+MeshSize{PointsOf{Volume{volMonopole[]};}} = (_use_fineMesh) ? lcMonopol : lcMonopol*4;
+MeshSize{PointsOf{Surface{surfaceAirInf[]};}} = (_use_fineMesh) ? lcair/2 : lcair;
+
+If(ang_section_degrees < 360)
+ Transfinite Line {9,12,24} = (_use_fineMesh) ? 21 : 11; // div length dipole
+ Transfinite Line {19,21,23,22} = (_use_fineMesh) ? 5 : 3; // div radius
+ Transfinite Line {6,11} = ang_section_degrees/45*((_use_fineMesh) ? 6 : 3); // lines cylinder
+ Transfinite Surface {4,5}; // cuts A & B
+ Transfinite Surface {11}; // cylindrical side
+ Transfinite Surface {2} = {12,6,5}; // top of terminal
+ Transfinite Surface {3} = {13,7,8}; // terminal
+ Transfinite Volume {volMonopole[]};
+EndIf
+
+// cut for the H formulation
+Cohomology(1) {{NONTERM},{}}; // this will generate the cut needed by the formulation in H, with no NONTERM + 1
diff --git a/ECE_full_wave_E-H_monopole/monopole.geo b/ECE_full_wave_E-H_monopole/monopole.geo
new file mode 100644
index 0000000000000000000000000000000000000000..f316d3e12b012e0fe6b5168d0d67e311455f0f57
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/monopole.geo
@@ -0,0 +1,69 @@
+Include "monopole_data.pro";
+
+If(Flag_3Dmodel==0) // 2D model
+ If (d == 0) // NO GAP
+ If (_use_pml == 0) // no PML
+ If (Flag_bndShape == 1) // cylindrical bnd
+ Printf("2.5D, cylindrical boundary, ECE radiant"); // 2.5D, cylindrical, ECE radiant
+ Include "./geo/monopole2D_cyl_rad.geo";
+ ElseIf (Flag_bndShape == 2) // spherical bnd
+ Printf("2.5D, spherical boundary, ECE radiant"); // 2.5D, spherical, ECE radiant
+ Include "./geo/monopole2D_sph_rad.geo";
+ Else // capsular bnd
+ Printf("2.5D, capsular boundary, ECE radiant"); // 2.5D, capsular, ECE radiant
+ Include "./geo/monopole2D_cap_rad.geo";
+ EndIf
+ Else // with PML
+ If (Flag_bndShape == 1) // cylindrical bnd
+ Printf("2.5D, cylindrical boundary, PML"); // 2.5D, cylindrical, PML
+ Include "./geo/monopole2D_cyl_pml.geo";
+ ElseIf (Flag_bndShape == 2) // spherical bnd
+ Printf("2D, spherical, PML");
+ Include "./geo/monopole2D_sph_pml.geo"; // 2.5D, spherical, PML
+ Else // capsular bnd
+ Printf("2.5D, capsular boundary, PML"); // 2.5D, capsular, PML
+ Include "./geo/monopole2D_cap_pml.geo";
+ EndIf
+ EndIf
+ Else
+ Printf("2.5D, spherical boundary, ECE radiant, gap > 0 or coax feeding"); // 2.5D, spherical, ECE radiant, gap or coax feeding
+ If ((Flag_bndShape == 2)&(_use_pml == 0))
+ If (d > 0) // gap (hard source feeding)
+ Include "./geo/monopole2Dgap_sph_rad.geo";
+ Else
+ Include "./geo/monopole2Dcoax_sph_rad.geo";
+ EndIf
+ Else
+ Printf("coax and gap, implemented in 2Daxi only with a circular boundary and radECE");
+ EndIf
+ EndIf
+Else // 3D model - only no gap, with spherical bnd
+ If (d == 0) // NO GAP
+ If (_use_pml == 0) // no PML
+ If (Flag_bndShape == 1) // cylindrical bnd
+ Printf("no gap, 3D, radiant ECE, cylindrical bnd - not implemented");
+ //Include "./geo/monopole3D_cyl_rad.geo";
+ ElseIf (Flag_bndShape == 2) // spherical bnd
+ Printf("no gap, 3D, radiant ECE, spherical bnd"); // 3D, spherical, ECE radiant
+ //Include "./geo/monopole3D_sph_rad.geo";
+ Include "./geo/monopole3D_sph_rad_new.geo";
+ Else // capsular bnd
+ Printf("no gap, 3D, radiant ECE, capsular bnd - not implemented");
+ //Include "./geo/monopole3D_cap_rad.geo";
+ EndIf
+ Else // with PML
+ If (Flag_bndShape == 1) // cylindrical bnd
+ Printf("no gap, 3D, PML, cylindrical bnd - not implemented");
+ //Include "./geo/monopole3D_cyl_pml.geo";
+ ElseIf (Flag_bndShape == 2) // spherical bnd
+ Printf("no gap, 3D, PML, spherical bnd");
+ Include "./geo/monopole3D_sph_pml.geo";
+ Else // capsular bnd
+ Printf("no gap, 3D, PML, capsular bnd - not implemented");
+ //Include "./geo/monopole3D_cap_pml.geo";
+ EndIf
+ EndIf
+ Else
+ Printf("coax and gap, not implemented in 3D");
+ EndIf
+EndIf
diff --git a/ECE_full_wave_E-H_monopole/monopole.pro b/ECE_full_wave_E-H_monopole/monopole.pro
new file mode 100644
index 0000000000000000000000000000000000000000..753c4da6fdcc9dd3d9058ff446487d53ffe0827d
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/monopole.pro
@@ -0,0 +1,11 @@
+Include "monopole_data.pro"
+
+If(Flag_3Dmodel==0) // 2D axi model
+ Include "monopole2Daxi.pro";
+ElseIf (d==0)
+ Include "monopole3D.pro";
+Else
+ Printf("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
+ Printf("!!!!!!!!!!!!!!!! coax and gap, not implemented in 3D !!!!!!!!!!!!!!!!");
+ Printf("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
+EndIf
diff --git a/ECE_full_wave_E-H_monopole/monopole2Daxi.pro b/ECE_full_wave_E-H_monopole/monopole2Daxi.pro
new file mode 100644
index 0000000000000000000000000000000000000000..f6bd266a6903b4a59828719ba8bde8f4e720afbc
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/monopole2Daxi.pro
@@ -0,0 +1,416 @@
+
+Group{
+ // Definition of the interface between the mesh file and GetDP
+ // Surfaces
+ Monopole = Region[{ MONOPOLE }] ;
+ Air = Region[{ AIR }] ;
+ Pml = Region[{}] ;
+ If (_use_pml)
+ Pml += Region[{ PML }] ;
+ EndIf
+ If (d < 0)
+ CoaxCond = Region[{ COAXCONDUCTOR }];
+ CoaxDiel = Region[{ COAXDIELECTRIC}];
+ EndIf
+
+ // Boundary
+ Terminal = Region[{ TERMINAL }] ;
+ Cut1H = Region[{ CUT }] ;
+
+ Axis = Region[{ AXIS }] ;
+ Ground = Region[{ GROUND }] ;
+
+ Sur_Terminals_FWece = Region[{Ground, Terminal}]; // all purely conducting terminals
+
+ // (Default: Flag_AnalysisTypeExtBnd==1) conducting ECE, fictitious boundary set to ground terminal - ~PEC
+ BoundaryNotTerminal = Region[{Cut1H}];
+ SurfAirInf = Region[{}];
+ If (Flag_AnalysisTypeExtBnd != 1)
+ SurfAirInf += Region[{ SURFAIRINF }] ;
+ BoundaryNotTerminal += Region[{SurfAirInf}];
+ EndIf
+
+ Sur_FW = Region[{Sur_Terminals_FWece, BoundaryNotTerminal}]; //all boundary
+ If(_use_radiantECE)
+ Sur_Terminals_FWeceAndInf = Region[{Sur_Terminals_FWece, SurfAirInf}];
+ EndIf
+
+ Vol_FW = Region[{Monopole,Air}];
+ If (d<0)
+ Vol_FW += Region[{CoaxCond,CoaxDiel}];
+ EndIf
+ If (_use_pml)
+ Vol_FW += Region[{Pml}];
+ EndIf
+ Dom_FW = Region[ {Vol_FW, Sur_FW, Axis} ];
+}
+
+Function{
+ // Material properties are defined piecewise, for elementary groups
+ // the values in the right hand side are defined in monopole_data.pro
+
+ // monopole
+ nu[#{Monopole}] = nuMonopole; // magnetic reluctivity [m/H]
+ epsilon[#{Monopole}] = epsMonopole; // permittivity [F/m]
+ sigma[#{Monopole}] = sigmaMonopole; // conductivity [S/m]
+ If (d>=0)
+ nu[#{Terminal}] = nuMonopole;
+ epsilon[#{Terminal}] = epsMonopole;
+ sigma[#{Terminal}] = sigmaMonopole;
+ EndIf
+
+ // air
+ nu[#{Air,Ground}] = nuAir;
+ epsilon[#{Air,Ground}] = epsAir;
+ sigma[#{Air,Ground}] = sigmaAir;
+ If (d>=0)
+ nu[#{Cut1H}] = nuAir;
+ epsilon[#{Cut1H}] = epsAir;
+ sigma[#{Cut1H}] = sigmaAir;
+ EndIf
+
+ // functions needed for PML only - begin
+ k[]= k0;
+ If (Flag_bndShape == 2) // SPHERICAL DOMAIN
+ Printf('=========================================> PML - (2) spherical boundary');
+ R[] = Sqrt[X[]^2 + Y[]^2 + Z[]^2];
+ rLoc[] = R[]-rb;
+ absFuncS[] = 1/(PmlDelta-rLoc[]);
+ absFuncF[] = -Log[1-rLoc[]/PmlDelta];
+
+ s1[] = Complex[1,-absFuncS[]/k[]] ; //cR
+ s2[] = Complex[1,-(1/R[])*absFuncF[]/k[]] ; //cStretch
+ nVec[] = XYZ[]/R[];
+ tVec[] = nVec[] /\ Vector[0,0,1];
+ nTen[] = SquDyadicProduct[nVec[]];
+ tTen[] = SquDyadicProduct[tVec[]];
+ pmlScal[] = s1[]*s2[];
+ pmlTens[] = (s2[]/s1[]) * nTen[] + (s1[]/s2[]) * tTen[];
+ tens[] = TensorSym[CompXX[pmlTens[]], CompXY[pmlTens[]], 0,
+ CompYY[pmlTens[]], 0,
+ (CompZZ[pmlTens[]]!=0) ? CompZZ[pmlTens[]]: pmlScal[] ];
+ ElseIf (Flag_bndShape == 1) // CYLINDRICAL DOMAIN
+ Printf('=========================================> PML - (1) Cylindrical boundary');
+ xLoc[] = Fabs[X[]]-xb;
+ yLoc[] = Fabs[Y[]]-yb;
+ zLoc[] = Fabs[Z[]]-zb;
+ DampingProfileX[] = (xLoc[]>=0) ? 1 / (PmlDelta-xLoc[]) : 0 ;
+ DampingProfileY[] = (yLoc[]>=0) ? 1 / (PmlDelta-yLoc[]) : 0 ;
+ DampingProfileZ[] = (zLoc[]>=0) ? 1 / (PmlDelta-zLoc[]) : 0 ;
+
+ cX[] = Complex[1,-DampingProfileX[]/k0] ;
+ cY[] = Complex[1,-DampingProfileY[]/k0] ;
+ cZ[] = (Flag_3Dmodel==1) ? Complex[1,-DampingProfileZ[]/k0] : 1. ;
+
+ t11[] = cY[]*cZ[]/cX[];
+ t22[] = cX[]*cZ[]/cY[];
+ t33[] = cX[]*cY[]/cZ[] ;
+ tens[] = TensorDiag[ t11[], t22[], t33[] ] ;
+ Else // CAPSULAR DOMAIN
+ Printf('=========================================> PML - (3) capsular boundary');
+ R[] = Sqrt[X[]^2 + (Y[]-L)^2 + Z[]^2];
+ rLoc[] = R[]-rb;
+ absFuncS[] = 1/(PmlDelta-rLoc[]);
+ absFuncF[] = -Log[1-rLoc[]/PmlDelta];
+
+ s1[] = Complex[1,-absFuncS[]/k[]] ; //cR
+ s2[] = Complex[1,-1/R[]*absFuncF[]/k[]] ; //cStretch
+
+ xLoc[] = X[]-rb;
+ absFuncSX[] = (xLoc[]>=0) ? 1/(PmlDelta-xLoc[]) : 0 ;
+ s1X[] = Complex[1,-absFuncSX[]/k[]] ; //cX[]
+
+ yLoc[] = Y[]-rb;
+ absFuncSY[] = (yLoc[]>=0) ? 1/(PmlDelta-yLoc[]) : 0 ;
+ s1Y[] = Complex[1,-absFuncSY[]/k[]] ; //cY[]
+
+ nVec[] = XYZ[]/R[];
+ tVec[] = nVec[] /\ Vector[0,0,1];
+ nTen[] = SquDyadicProduct[nVec[]];
+ tTen[] = SquDyadicProduct[tVec[]];
+ pmlScal[] = s1[]*s2[];
+ pmlTens[] = (s2[]/s1[]) * nTen[] + (s1[]/s2[]) * tTen[];
+
+ tens[] = TensorSym[(Fabs[Y[]]>=L) ? CompXX[pmlTens[]]: s1Y[]/s1X[],
+ (Fabs[Y[]]>=L) ? CompXY[pmlTens[]]: 0,
+ 0,
+ (Fabs[Y[]]>=L) ? CompYY[pmlTens[]]: s1X[]/s1Y[],
+ 0,
+ (Fabs[Y[]]>=L) ? ((CompZZ[pmlTens[]]!=0) ? CompZZ[pmlTens[]]:pmlScal[]) : s1X[]*s1Y[]];
+ EndIf
+
+ epsilon[ #{Pml} ] = eps0 * tens[] ;
+ nu[ #{Pml} ] = nu0 / tens[] ;
+
+ epsilon[ #{SurfAirInf} ] = eps0 ;
+ nu[ #{SurfAirInf} ] = nu0;
+
+ t11d[] = 1;
+ sigma[ #{Pml,SurfAirInf} ] = sigmaAir * TensorDiag[ t11d[], t11d[], t11d[] ];
+
+ // coax
+ If (d < 0)
+ nu[#{CoaxCond}] = nuCoaxCond;
+ epsilon[#{CoaxCond}] = epsCoaxCond;
+ sigma[#{CoaxCond}] = sigmaCoaxCond;
+
+ nu[#{CoaxDiel}] = nuCoaxDiel;
+ epsilon[#{CoaxDiel}] = epsCoaxDiel;
+ sigma[#{CoaxDiel}] = sigmaCoaxDiel;
+ EndIf
+
+ mu[] = 1/nu[];
+
+ epsilonNu[] = epsilon[]*nu[] ;
+
+ // excitation - voltage or current
+ // No need of Flag_AnalysisType here => constraints
+ VTerminal1[] = (Flag_AnalysisTypeFormulation==0) ? -1 : 1;
+ ITerminal1[] = (Flag_AnalysisTypeFormulation==0) ? -1 : 1;
+}
+
+ Constraint{
+
+ // voltage excited terminals, formulation in E
+ { Name SetTerminalPotential; Type Assign;
+ Case {
+ { Region Ground; Value 0.; }
+ If(Flag_AnalysisType == 0)
+ { Region Terminal; Value VTerminal1[]; }
+ EndIf
+ }
+ }
+
+ // current excited terminals, formulation in E
+ { Name SetTerminalCurrent; Type Assign;
+ Case {
+ If(Flag_AnalysisType == 1)
+ { Region Terminal; Value ITerminal1[]/h2Ddepth; }
+ EndIf
+ }
+ }
+
+ // current excited terminals, formulation in H, constraints are on the cut
+ { Name SetTerminalCurrentH; Type Assign;
+ Case {
+ If(Flag_AnalysisType == 1)
+ { Region Cut1H; Value ITerminal1[]/h2Ddepth; }
+ EndIf
+ }
+ }
+
+ // voltage excited terminals, formulation in H, constraints are on the cut
+ { Name SetTerminalPotentialH ;
+ Case {
+ If(Flag_AnalysisType == 0)
+ { Region Cut1H; Value VTerminal1[] ; }
+ EndIf
+ }
+ }
+
+ // formulation in H, 2.5D axi, magnetic field strength has to be imposed as zero on the axis
+ { Name ImposeHonAxis ;
+ Case {
+ If(Flag_AnalysisTypeFormulation == 1 && Flag_Axi == 1) // H
+ { Region Axis; Value 0 ; }
+ EndIf
+ }
+ }
+
+ // formulation in E, classical PEC on the fictitious boundary (classical, used only for some checks)
+ { Name ElectricVoltages ;
+ Case {
+ If(Flag_AnalysisTypeFormulation == 0) // E
+ If (_use_radiantECE) // here it means no pml
+ If (Flag_SilverMuller == 0) // no radiant
+ If (Flag_PEC)
+ { Region SurfAirInf; Value 0 ; } //Et will be imposed to zero on the fictitious boundary
+ EndIf
+ EndIf
+ EndIf
+ EndIf
+ }
+ }
+
+ // formulation in H, classical PMC on the fictitious boundary (classical, used only for some checks
+ { Name MagneticVoltages ;
+ Case {
+ If(Flag_AnalysisTypeFormulation==1) // H
+ If (_use_radiantECE)
+ If (Flag_SilverMuller == 0)
+ If (Flag_PEC==0) // this means PMC
+ { Region SurfAirInf; Value 0 ; } // radiant, classical, PMC
+ EndIf
+ EndIf
+ EndIf
+ EndIf
+ }
+ }
+}
+
+modelDir = "../";
+// folders where the frequency responses ins Touchstone files will be saved
+If (Flag_AnalysisType == 0) // voltage excitation
+ If (Flag_AnalysisTypeFormulation == 0)
+ Dir = "res/FWeceBC_voltExc/"; // E
+ Else
+ Dir = "res/H-FWeceBC_voltExc/"; // H
+ EndIf
+Else // current excitation
+ If (Flag_AnalysisTypeFormulation ==0)
+ Dir = "res/FWeceBC_crtExc/"; // E
+ Else
+ Dir = "res/H-FWeceBC_crtExc/"; // H
+ EndIf
+EndIf
+
+// The formulation and its tools
+If (Flag_AnalysisTypeFormulation == 0) // E
+ Include "./pro_problemIndependent/FullWave_E_eceRadiant_SISO_secondOrder.pro"
+Else // H, 2.5D
+ Include "./pro_problemIndependent/FullWave_H_eceRadiant_SISO_cplx_2D_and_AXI_secondOrder.pro"
+EndIf
+
+/* Output results */
+Macro Change_post_options
+Echo[Str[ "nv = PostProcessing.NbViews;",
+ "For v In {0:nv-1}",
+ "View[v].NbIso = 25;",
+ "View[v].RangeType = 3;" ,// per timestep
+ "View[v].ShowTime = 3;",
+ "View[v].IntervalsType = 3;",
+ "EndFor"
+ ], File "post.opt"];
+Return
+
+Macro WriteZ_header // useful for current excitation
+ If(Flag_AnalysisTypeFormulation == 0) // E
+ Echo[ Str[Sprintf("# Hz Z RI R %g", 50)], Format Table, File > StrCat[Dir, "test_Z_RI_formE.s1p"] ];
+ Else // H
+ Echo[ Str[Sprintf("# Hz Z RI R %g", 50)], Format Table, File > StrCat[Dir, "test_Z_RI_formH.s1p"] ];
+ EndIf
+Return
+
+Macro WriteY_header // useful for voltage excitation
+ If(Flag_AnalysisTypeFormulation == 0) // E
+ Echo[ Str[Sprintf("# Hz Y RI R %g", 50)], Format Table, File > StrCat[Dir, "test_Y_RI_formE.s1p"] ];
+ Else // H
+ Echo[ Str[Sprintf("# Hz Y RI R %g", 50)], Format Table, File > StrCat[Dir, "test_Y_RI_formH.s1p"] ];
+ EndIf
+Return
+
+PostOperation {
+ If (Flag_AnalysisTypeFormulation == 0) // formulation in E
+
+ { Name TransferMatrix; NameOfPostProcessing FW_E_ece ;
+ Operation {
+ If(Flag_AnalysisType == 1) // current excitation
+ If (Freq == freqs(0)) // write the headear in the s1p file
+ Call WriteZ_header;
+ Printf("======================> Writting header with first Freq=%g",Freq);
+ EndIf
+ Print [ Vterminals, OnRegion Terminal, Format Table , File > StrCat[Dir, "test_Z_RI_formE.s1p"]] ;
+ Else // voltage excitation
+ If (Freq == freqs(0)) // write the headear in the s1p file
+ Call WriteY_header;
+ Printf("======================> Writting header with first Freq=%g",Freq);
+ EndIf
+ Print [ I, OnRegion Terminal, Format Table , File > StrCat[Dir, "test_Y_RI_formE.s1p"]] ;
+ EndIf
+
+ If(Flag_AnalysisType==1) // current excitation
+ Print[ Rin, OnRegion Terminal, Format Table, SendToServer StrCat[pp0,"R=re(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempRin.dat"] ] ;
+ Print[ Xin, OnRegion Terminal, Format Table, SendToServer StrCat[pp1,"X=im(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempXin.dat"] ] ;
+ Print[ AbsZin, OnRegion Terminal, Format Table, SendToServer StrCat[pp2,"|Z|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsZin.dat"] ];
+ Print[ ArgZin, OnRegion Terminal, Format Table, SendToServer StrCat[pp3,"arg(Z) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgZin.dat"] ] ;
+ Else // voltage excitation
+ Print[ Yin, OnRegion Terminal, Format Table, SendToServer StrCat[pp0,"G=re(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempGin.dat"] ] ;
+ Print[ Yin, OnRegion Terminal, Format Table, SendToServer StrCat[pp1,"B=im(Y)"]{1}, Color "Ivory", File StrCat[Dir,"tempBin.dat"] ] ;
+
+ Print[ AbsYin, OnRegion Terminal, Format Table, SendToServer StrCat[pp2,"|Y|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsYin.dat"] ];
+ Print[ ArgYin, OnRegion Terminal, Format Table, SendToServer StrCat[pp3,"arg(Y) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgYin.dat"] ] ;
+ EndIf
+ }
+ }
+
+ { Name Maps; NameOfPostProcessing FW_E_ece ;
+ Operation {
+ Print [ gradV, OnElementsOf ElementsOf[Vol_FW,OnOneSideOf Sur_FW], File StrCat[Dir,"map_gradV.pos" ]];
+ Print [ Vterminals, OnElementsOf ElementsOf[Sur_FW], File StrCat[Dir,"map_Vterminals.pos" ]];
+
+ Print [ E, OnElementsOf Vol_FW, File StrCat[Dir,"map_E.pos" ]];
+ Print [ J, OnElementsOf Vol_FW, File StrCat[Dir,"map_J.pos" ]];
+ Print [ rmsJ, OnElementsOf Vol_FW, File StrCat[Dir,"map_absJ.pos" ]];
+
+ Print [ H, OnElementsOf Region[{Vol_FW,-Pml}] , File StrCat[Dir,"map_H.pos" ]];
+
+ Print [ VnotTerminals, OnElementsOf Vol_FW, File StrCat[Dir,"map_VnotTerminals.pos" ]];
+
+ If(Flag_AnalysisType == 1) // current excitation
+ Print[ Rin, OnRegion Terminal, Format Table, SendToServer StrCat[pp0,"R=re(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempRin.dat"] ] ;
+ Print[ Xin, OnRegion Terminal, Format Table, SendToServer StrCat[pp1,"X=im(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempXin.dat"] ] ;
+ Print[ AbsZin, OnRegion Terminal, Format Table, SendToServer StrCat[pp2,"|Z|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsZin.dat"] ];
+ Print[ ArgZin, OnRegion Terminal, Format Table, SendToServer StrCat[pp3,"arg(Z) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgZin.dat"] ] ;
+ Else // voltage excitation
+ Print[ Yin, OnRegion Terminal, Format Table, SendToServer StrCat[pp0,"G=re(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempGin.dat"] ] ;
+ Print[ Yin, OnRegion Terminal, Format Table, SendToServer StrCat[pp1,"B=im(Y)"]{1}, Color "Ivory", File StrCat[Dir,"tempBin.dat"] ] ;
+
+ Print[ AbsYin, OnRegion Terminal, Format Table, SendToServer StrCat[pp2,"|Y|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsYin.dat"] ];
+ Print[ ArgYin, OnRegion Terminal, Format Table, SendToServer StrCat[pp3,"arg(Y) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgYin.dat"] ] ;
+ EndIf
+ Call Change_post_options;
+ }
+ }
+ Else // formulation in H
+ { Name TransferMatrix; NameOfPostProcessing FW_H_ece ;
+ Operation {
+ If(Flag_AnalysisType == 1) // current excitation
+ If (Freq == freqs(0)) // write the headear in the s1p file
+ Call WriteZ_header;
+ Printf("======================> Writting header with first Freq=%g",Freq);
+ EndIf
+ Print [ Vterminals, OnRegion Cut1H, Format Table , File > StrCat[Dir, "test_Z_RI_formH.s1p"]] ;
+ Else // voltage excitation
+ If (Freq == freqs(0)) // write the headear in the s1p file
+ Call WriteY_header;
+ Printf("======================> Writting header with first Freq=%g",Freq);
+ EndIf
+ Print [ I, OnRegion Cut1H, Format Table , File > StrCat[Dir, "test_Y_RI_formH.s1p"]] ;
+ EndIf
+
+ If(Flag_AnalysisType == 1) // current excitation
+ Print[ Rin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp0,"R=re(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempRin.dat"] ] ;
+ Print[ Xin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp1,"X=im(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempXin.dat"] ] ;
+ Print[ AbsZin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp2,"|Z|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsZin.dat"] ];
+ Print[ ArgZin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp3,"arg(Z) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgZin.dat"] ] ;
+ Else // voltage excitation
+ Print[ Gin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp0,"G=re(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempGin.dat"] ] ;
+ Print[ Bin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp1,"B=im(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempBin.dat"] ] ;
+ Print[ AbsYin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp2,"|Y|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsYin.dat"] ] ;
+ Print[ ArgYin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp3,"arg(Y) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgYin.dat"] ] ;
+ EndIf
+ }
+ }
+
+ { Name Maps; NameOfPostProcessing FW_H_ece ;
+ Operation {
+ Print [ E, OnElementsOf Region[{Vol_FW,-Pml}], File StrCat[Dir,"map_E.pos" ]];
+ Print [ J, OnElementsOf Region[{Vol_FW,-Pml}], File StrCat[Dir,"map_J.pos" ]];
+ Print [ rmsJ, OnElementsOf Region[{Vol_FW,-Pml}], File StrCat[Dir,"map_absJ.pos" ]];
+
+ Print [ H, OnElementsOf Vol_FW, File StrCat[Dir,"map_H.pos" ]];
+
+ If(Flag_AnalysisType == 1) // current excitation
+ Print[ Rin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp0,"R=re(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempRin.dat"] ] ;
+ Print[ Xin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp1,"X=im(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempXin.dat"] ] ;
+ Else // voltage excitation
+ Print[ Gin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp0,"G=re(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempGin.dat"] ] ;
+ Print[ Bin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp1,"B=im(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempBin.dat"] ] ;
+ Print[ AbsYin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp2,"|Y|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsYin.dat"] ] ;
+ Print[ ArgYin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp3,"arg(Y) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgYin.dat"] ] ;
+ EndIf
+ Call Change_post_options;
+ }
+ }
+ EndIf
+}
diff --git a/ECE_full_wave_E-H_monopole/monopole3D.pro b/ECE_full_wave_E-H_monopole/monopole3D.pro
new file mode 100644
index 0000000000000000000000000000000000000000..0f4bc3338f12b95ac8aee2ee626305d7e8df4309
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/monopole3D.pro
@@ -0,0 +1,431 @@
+Group{
+
+ // Surfaces
+ Monopole = Region[{ MONOPOLE }] ;
+ Air = Region[{ AIR }] ;
+ Pml = Region[{}] ;
+ If (_use_pml)
+ Pml += Region[{ PML }] ;
+ EndIf
+
+ If (d < 0)
+ CoaxCond = Region[{ COAXCONDUCTOR }];
+ CoaxDiel = Region[{ COAXDIELECTRIC}];
+ EndIf
+
+ // Boundary
+ Terminal = Region[{ TERMINAL }] ;
+ NonTerm = Region[{ NONTERM }];
+
+ Cut1H = Region[{}] ;
+ If(Flag_AnalysisTypeFormulation == 1) // H
+ Cut1H += Region[{ CUT1H }] ;
+ EndIf
+
+ Ground = Region[{ GROUND }] ;
+
+ Sur_Terminals_FWece = Region[{Ground, Terminal}]; // all purely conducting terminals
+
+ // (Default: Flag_AnalysisTypeExtBnd==1) conducting ECE, fictitious boundary set to ground terminal - ~PEC
+ BoundaryNotTerminal = Region[{NonTerm}];
+ SurfAirInf = Region[{}];
+ If (Flag_AnalysisTypeExtBnd != 1)
+ SurfAirInf += Region[{ SURFAIRINF }] ;
+ BoundaryNotTerminal += Region[{SurfAirInf}];
+ EndIf
+
+ Sur_FW = Region[{Sur_Terminals_FWece, BoundaryNotTerminal}]; //all boundary
+ If(_use_radiantECE)
+ Sur_Terminals_FWeceAndInf = Region[{Sur_Terminals_FWece, SurfAirInf}];
+ EndIf
+
+ Vol_FW = Region[{Monopole,Air}];
+ If (d<0)
+ Vol_FW += Region[{CoaxCond,CoaxDiel}];
+ EndIf
+
+ If (_use_pml)
+ Vol_FW += Region[{Pml}];
+ EndIf
+ Dom_FW = Region[ {Vol_FW, Sur_FW} ];
+}
+
+Function{
+ // Material properties are defined piecewise, for elementary groups
+ // the values in the right hand side are defined in monopole_data.pro
+
+ // monopole
+ nu[#{Monopole}] = nuMonopole; // magnetic reluctivity [m/H]
+ epsilon[#{Monopole}] = epsMonopole; // permittivity [F/m]
+ sigma[#{Monopole}] = sigmaMonopole; // conductivity [S/m]
+ If (d>=0)
+ nu[#{Terminal}] = nuMonopole;
+ epsilon[#{Terminal}] = epsMonopole;
+ sigma[#{Terminal}] = sigmaMonopole;
+ EndIf
+
+ // air
+ nu[#{Air,Ground}] = nuAir;
+ epsilon[#{Air,Ground}] = epsAir;
+ sigma[#{Air,Ground}] = sigmaAir;
+ If (d>=0)
+ nu[#{NonTerm}] = nuAir;
+ epsilon[#{NonTerm}] = epsAir;
+ sigma[#{NonTerm}] = sigmaAir;
+ EndIf
+
+ // functions needed for PML only - begin
+ k[]= k0;
+ If (Flag_bndShape == 2) // SPHERICAL DOMAIN
+ Printf('=========================================> PML - (2) spherical boundary');
+ R[] = Sqrt[X[]^2 + Y[]^2 + Z[]^2];
+ rLoc[] = R[]-rb;
+ absFuncS[] = (rLoc[]>=0) ? 1/(PmlDelta-rLoc[]) : 0.;
+ absFuncF[] = -Log[1-rLoc[]/PmlDelta];
+
+ s1[] = Complex[1,-absFuncS[]/k[]] ;
+ s2[] = Complex[1,-1/R[]*absFuncF[]/k[]] ;
+ nVec[] = XYZ[]/R[];
+ tVec[] = nVec[] /\ Vector[0,0,1];
+ nTen[] = SquDyadicProduct[nVec[]];
+ tTen[] = SquDyadicProduct[tVec[]];
+ pmlScal[] = s1[]*s2[];
+ pmlTens[] = (s2[]/s1[]) * nTen[] + (s1[]/s2[]) * tTen[];
+ dettTen[] = Det[tTen[]];
+
+ tens[] = (dettTen[]==0) ?
+ TensorSym[(CompXX[tTen[]]!=0) ? CompXX[pmlTens[]]: pmlScal[], CompXY[pmlTens[]], CompXZ[pmlTens[]],
+ (CompYY[tTen[]]!=0) ? CompYY[pmlTens[]]: pmlScal[], CompYZ[pmlTens[]],
+ (CompZZ[tTen[]]!=0) ? CompZZ[pmlTens[]]: pmlScal[] ] : pmlTens[] ;
+
+ ElseIf (Flag_bndShape == 1) // CYLINDRICAL DOMAIN
+ Printf('=========================================> (1) Cylindrical boundary');
+ xLoc[] = Fabs[X[]]-xb;
+ yLoc[] = Fabs[Y[]]-yb;
+ zLoc[] = Fabs[Z[]]-zb;
+ DampingProfileX[] = (xLoc[]>=0) ? 1 / (PmlDelta-xLoc[]) : 0 ;
+ DampingProfileY[] = (yLoc[]>=0) ? 1 / (PmlDelta-yLoc[]) : 0 ;
+ DampingProfileZ[] = (zLoc[]>=0) ? 1 / (PmlDelta-zLoc[]) : 0 ;
+
+ cX[] = Complex[1,-DampingProfileX[]/k0] ;
+ cY[] = Complex[1,-DampingProfileY[]/k0] ;
+ cZ[] = Complex[1,-DampingProfileZ[]/k0] ;
+
+ t11[] = cY[]*cZ[]/cX[];
+ t22[] = cX[]*cZ[]/cY[];
+ t33[] = cX[]*cY[]/cZ[] ;
+ tens[] = TensorDiag[ t11[], t22[], t33[] ] ;
+
+ Else // CAPSULAR DOMAIN
+ Printf('=========================================> (3) capsular boundary');
+ Ldipole = 2*L; // L is the monopole length so that not to change the code below (for the momment)
+ R[] = Sqrt[X[]^2 + Y[]^2 + Z[]^2]; // No need of moving the origin as we used the actual normal and tangent
+ rLoc[] = R[]-rb;
+ absFuncS[] = 1/(PmlDelta-rLoc[]);
+ absFuncF[] = -Log[1-rLoc[]/PmlDelta];
+
+ s1[] = Complex[1,-absFuncS[]/k[]] ; //cR
+ s2[] = Complex[1,-1/R[]*absFuncF[]/k[]] ; //cStretch
+
+ xLoc[] = X[]-rb;
+ absFuncSX[] = (xLoc[]>=0) ? 1/(PmlDelta-xLoc[]) : 0 ;
+ s1X[] = Complex[1,-absFuncSX[]/k[]] ; //cX[]
+
+ yLoc[] = Y[]-rb;
+ absFuncSY[] = (yLoc[]>=0) ? 1/(PmlDelta-yLoc[]) : 0 ;
+ s1Y[] = Complex[1,-absFuncSY[]/k[]] ; //cY[]
+
+ zLoc[] = Z[]-rb;
+ absFuncSZ[] = (zLoc[]>=0) ? 1/(PmlDelta-zLoc[]) : 0 ;
+ s1Z[] = Complex[1,-absFuncSZ[]/k[]] ; //cY[]
+
+ nVec[] = XYZ[]/R[];
+ tVec[] = (nVec[] /\ Vector[0,0,1]);
+ nTen[] = SquDyadicProduct[nVec[]];
+ tTen[] = SquDyadicProduct[tVec[]];
+ pmlScal[] = s1[]*s2[];
+ pmlTens[] = (s2[]/s1[]) * nTen[] + (s1[]/s2[]) * tTen[];
+
+ tens[] = TensorSym[ // ++++ TO CHECK
+ (Fabs[Y[]]>=Ldipole/2) ? CompXX[pmlTens[]]: s1Y[]*s1Z[]/s1X[],
+ (Fabs[Y[]]>=Ldipole/2) ? CompXY[pmlTens[]]: 0,
+ (Fabs[Y[]]>=Ldipole/2) ? CompXZ[pmlTens[]]: 0,
+ 0,
+ 0,
+ (Fabs[Y[]]>=Ldipole/2) ? ((CompZZ[pmlTens[]]!=0) ? CompZZ[pmlTens[]]:pmlScal[]) : s1X[]*s1Y[]/s1Z[]];
+ EndIf
+
+ epsilon[ #{Pml} ] = eps0 * tens[] ;
+ nu[ #{Pml} ] = nu0 / tens[] ;
+
+ epsilon[ #{SurfAirInf} ] = eps0 ;
+ nu[ #{SurfAirInf} ] = nu0 ;
+
+ t11d[] = 1;
+ sigma[ #{Pml,SurfAirInf} ] = sigmaAir * TensorDiag[ t11d[], t11d[], t11d[] ];
+
+ // coax
+ If (d < 0)
+ nu[#{CoaxCond,Terminal}] = nuCoaxCond;
+ epsilon[#{CoaxCond,Terminal}] = epsCoaxCond;
+ sigma[#{CoaxCond,Terminal}] = sigmaCoaxCond;
+
+ nu[#{CoaxDiel,NonTerm,GroundCoax}] = nuCoaxDiel;
+ epsilon[#{CoaxDiel,NonTerm,GroundCoax}] = epsCoaxDiel;
+ sigma[#{CoaxDiel,NonTerm,GroundCoax}] = sigmaCoaxDiel;
+ EndIf
+
+ mu[] = 1/nu[];
+ epsilonNu[] = epsilon[]*nu[];
+
+ // excitation - voltage or current
+ VTerminal1[] = (Flag_AnalysisTypeFormulation==0) ? -1 : 1;
+ ITerminal1[] = (Flag_AnalysisTypeFormulation==0) ? -1 : 1;
+}
+
+Constraint{
+
+ // voltage excited terminals, formulation in E
+ { Name SetTerminalPotential; Type Assign;
+ Case {
+ { Region Ground; Value 0.; }
+ If(Flag_AnalysisType == 0)
+ { Region Terminal; Value VTerminal1[]; }
+ EndIf
+ }
+ }
+
+ // current excited terminals, formulation in E
+ { Name SetTerminalCurrent; Type Assign;
+ Case {
+ If(Flag_AnalysisType == 1)
+ { Region Terminal; Value ITerminal1[]/h2Ddepth; }
+ EndIf
+ }
+ }
+
+ // current excited terminals, formulation in H, constraints are on the cut
+ { Name SetTerminalCurrentH; Type Assign;
+ Case {
+ If(Flag_AnalysisType == 1)
+ { Region Cut1H; Value ITerminal1[]/h2Ddepth; }
+ EndIf
+ }
+ }
+
+ // voltage excited terminals, formulation in H, constraints are on the cut
+ { Name SetTerminalPotentialH ;
+ Case {
+ If(Flag_AnalysisType == 0)
+ { Region Cut1H; Value VTerminal1[] ; }
+ EndIf
+ }
+ }
+
+ // formulation in H, 2.5D axi, magnetic field strength has to be imposed as zero on the axis
+ { Name ImposeHonAxis ;
+ Case {
+ If(Flag_AnalysisTypeFormulation == 1) // H
+ If(Flag_Axi == 1) // axi
+ { Region Axis; Value 0 ; }
+ EndIf
+ EndIf
+ }
+ }
+
+ // formulation in E, classical PEC on the fictitious boundary (classical, used only for some checks)
+ { Name ElectricVoltages ;
+ Case {
+ If(Flag_AnalysisTypeFormulation == 0) // E
+ If (_use_radiantECE) // here it means no pml
+ If (Flag_SilverMuller == 0) // no radiant
+ If (Flag_PEC)
+ { Region SurfAirInf; Value 0 ; } //Et will be imposed to zero on the fictitious boundary
+ EndIf
+ EndIf
+ EndIf
+ EndIf
+ }
+ }
+
+ // formulation in H, classical PMC on the fictitious boundary (classical, used only for some checks
+ { Name MagneticVoltages ;
+ Case {
+ If(Flag_AnalysisTypeFormulation==1) // H
+ If (_use_radiantECE)
+ If (Flag_SilverMuller == 0)
+ If (Flag_PEC==0) // this means PMC
+ { Region SurfAirInf; Value 0 ; } // radiant, classical, PMC
+ EndIf
+ EndIf
+ EndIf
+ EndIf
+ }
+ }
+}
+
+modelDir = "../";
+// folders where the frequency responses ins Touchstone files will be saved
+If (Flag_AnalysisType == 0) // voltage excitation
+ If (Flag_AnalysisTypeFormulation == 0)
+ Dir = "res/FWeceBC_voltExc/"; // E
+ Else
+ Dir = "res/H-FWeceBC_voltExc/"; // H
+ EndIf
+Else // current excitation
+ If (Flag_AnalysisTypeFormulation ==0)
+ Dir = "res/FWeceBC_crtExc/"; // E
+ Else
+ Dir = "res/H-FWeceBC_crtExc/"; // H
+ EndIf
+EndIf
+
+// The formulation and its tools
+ If (Flag_AnalysisTypeFormulation == 0) // E
+ Include "./pro_problemIndependent/FullWave_E_eceRadiant_SISO_secondOrder.pro"
+ Else // H, 2.5D
+ Include "./pro_problemIndependent/FullWave_H_eceRadiant_SISO_cplx_3D_secondOrder.pro"
+ EndIf
+
+// Output results
+Macro Change_post_options
+Echo[Str[ "nv = PostProcessing.NbViews;",
+ "For v In {0:nv-1}",
+ "View[v].NbIso = 25;",
+ "View[v].RangeType = 3;" ,// per timestep
+ "View[v].ShowTime = 3;",
+ "View[v].IntervalsType = 3;",
+ "EndFor"
+ ], File "post.opt"];
+Return
+
+Macro WriteZ_header // useful for current excitation
+ If(Flag_AnalysisTypeFormulation == 0) // E
+ Echo[ Str[Sprintf("# Hz Z RI R %g", 50)], Format Table, File > StrCat[Dir, "test_Z_RI_formE.s1p"] ];
+ Else // H
+ Echo[ Str[Sprintf("# Hz Z RI R %g", 50)], Format Table, File > StrCat[Dir, "test_Z_RI_formH.s1p"] ];
+ EndIf
+Return
+Macro WriteY_header // useful for voltage excitation
+ If(Flag_AnalysisTypeFormulation == 0) // E
+ Echo[ Str[Sprintf("# Hz Y RI R %g", 50)], Format Table, File > StrCat[Dir, "test_Y_RI_formE.s1p"] ];
+ Else // H
+ Echo[ Str[Sprintf("# Hz Y RI R %g", 50)], Format Table, File > StrCat[Dir, "test_Y_RI_formH.s1p"] ];
+ EndIf
+Return
+
+PostOperation {
+ If (Flag_AnalysisTypeFormulation == 0) // formulation in E
+
+ { Name TransferMatrix; NameOfPostProcessing FW_E_ece ;
+ Operation {
+ If(Flag_AnalysisType == 1) // current excitation
+ If (Freq == freqs(0)) // write the headear in the s1p file
+ Call WriteZ_header;
+ Printf("======================> Writting header with first Freq=%g",Freq);
+ EndIf
+ Print [ Vterminals, OnRegion Terminal, Format Table , File > StrCat[Dir, "test_Z_RI_formE.s1p"]] ;
+ Else // voltage excitation
+ If (Freq == freqs(0)) // write the headear in the s1p file
+ Call WriteY_header;
+ Printf("======================> Writting header with first Freq=%g",Freq);
+ EndIf
+ Print [ I, OnRegion Terminal, Format Table , File > StrCat[Dir, "test_Y_RI_formE.s1p"]] ;
+ EndIf
+
+ If(Flag_AnalysisType==1) // current excitation
+ Print[ Rin, OnRegion Terminal, Format Table, SendToServer StrCat[pp0,"R=re(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempRin.dat"] ] ;
+ Print[ Xin, OnRegion Terminal, Format Table, SendToServer StrCat[pp1,"X=im(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempXin.dat"] ] ;
+ Print[ AbsZin, OnRegion Terminal, Format Table, SendToServer StrCat[pp2,"|Z|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsZin.dat"] ];
+ Print[ ArgZin, OnRegion Terminal, Format Table, SendToServer StrCat[pp3,"arg(Z) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgZin.dat"] ] ;
+ Else // voltage excitation
+ Print[ Gin, OnRegion Terminal, Format Table, SendToServer StrCat[pp0,"G=re(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempGin.dat"] ] ;
+ Print[ Bin, OnRegion Terminal, Format Table, SendToServer StrCat[pp1,"B=im(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempBin.dat"] ] ;
+ Print[ AbsYin, OnRegion Terminal, Format Table, SendToServer StrCat[pp2,"|Y|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsYin.dat"] ];
+ Print[ ArgYin, OnRegion Terminal, Format Table, SendToServer StrCat[pp3,"arg(Y) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgYin.dat"] ] ;
+ EndIf
+ }
+ }
+
+ { Name Maps; NameOfPostProcessing FW_E_ece ;
+ Operation {
+ Print [ gradV, OnElementsOf ElementsOf[Vol_FW,OnOneSideOf Sur_FW], File StrCat[Dir,"map_gradV.pos" ]];
+ Print [ Vterminals, OnElementsOf Sur_FW, File StrCat[Dir,"map_Vterminals.pos" ]];
+
+ Print [ E, OnElementsOf Vol_FW, File StrCat[Dir,"map_E.pos" ]];
+ Print [ J, OnElementsOf Vol_FW, File StrCat[Dir,"map_J.pos" ]];
+ Print [ rmsJ, OnElementsOf Vol_FW, File StrCat[Dir,"map_absJ.pos" ]];
+ Print [ H, OnElementsOf Region[{Vol_FW,-Pml}] , File StrCat[Dir,"map_H.pos" ]];
+ Print [ VnotTerminals, OnElementsOf ElementsOf[Vol_FW, OnOneSideOf NonTerm], File StrCat[Dir,"map_VnotTerminals.pos" ]];
+
+ If(Flag_AnalysisType == 1) // current excitation
+ Print[ Rin, OnRegion Terminal, Format Table, SendToServer StrCat[pp0,"R=re(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempRin.dat"] ] ;
+ Print[ Xin, OnRegion Terminal, Format Table, SendToServer StrCat[pp1,"X=im(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempXin.dat"] ] ;
+ Print[ AbsZin, OnRegion Terminal, Format Table, SendToServer StrCat[pp2,"|Z|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsZin.dat"] ];
+ Print[ ArgZin, OnRegion Terminal, Format Table, SendToServer StrCat[pp3,"arg(Z) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgZin.dat"] ] ;
+ Else // voltage excitation
+ Print[ Gin, OnRegion Terminal, Format Table, SendToServer StrCat[pp0,"G=re(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempGin.dat"] ] ;
+ Print[ Bin, OnRegion Terminal, Format Table, SendToServer StrCat[pp1,"B=im(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempBin.dat"] ] ;
+ Print[ AbsYin, OnRegion Terminal, Format Table, SendToServer StrCat[pp2,"|Y|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsYin.dat"] ];
+ Print[ ArgYin, OnRegion Terminal, Format Table, SendToServer StrCat[pp3,"arg(Y) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgYin.dat"] ] ;
+ EndIf
+ Call Change_post_options;
+ }
+ }
+
+ Else // formulation in H
+
+ { Name TransferMatrix; NameOfPostProcessing FW_H_ece ;
+ Operation {
+ If(Flag_AnalysisType == 1) // current excitation
+ If (Freq == freqs(0)) // write the headear in the s1p file
+ Call WriteZ_header;
+ Printf("======================> Writting header with first Freq=%g",Freq);
+ EndIf
+ Print [ Vterminals, OnRegion Cut1H, Format Table , File > StrCat[Dir, "test_Z_RI_formH.s1p"]] ;
+ Else // voltage excitation */
+ If (Freq==freqs(0)) // write the headear in the s1p file
+ Call WriteY_header;
+ Printf("======================> Writting header with first Freq=%g",Freq);
+ EndIf
+ Print [ I, OnRegion Cut1H, Format Table , File > StrCat[Dir, "test_Y_RI_formH.s1p"]] ;
+ EndIf
+
+ If(Flag_AnalysisType == 1) // current excitation
+ Print[ Rin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp0,"R=re(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempRin.dat"] ] ;
+ Print[ Xin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp1,"X=im(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempXin.dat"] ] ;
+ Print[ AbsZin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp2,"|Z|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsZin.dat"] ];
+ Print[ ArgZin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp3,"arg(Z) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgZin.dat"] ] ;
+ Else // voltage excitation
+ Print[ Gin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp0,"G=re(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempGin.dat"] ] ;
+ Print[ Bin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp1,"B=im(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempBin.dat"] ] ;
+ Print[ AbsYin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp2,"|Y|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsYin.dat"] ] ;
+ Print[ ArgYin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp3,"arg(Y) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgYin.dat"] ] ;
+ EndIf
+ }
+ }
+
+ { Name Maps; NameOfPostProcessing FW_H_ece ;
+ Operation {
+ Print [ E, OnElementsOf Region[{Vol_FW,-Pml}], File StrCat[Dir,"map_E.pos" ]];
+ Print [ J, OnElementsOf Region[{Vol_FW,-Pml}], File StrCat[Dir,"map_J.pos" ]];
+ Print [ rmsJ, OnElementsOf Region[{Vol_FW,-Pml}], File StrCat[Dir,"map_absJ.pos" ]];
+
+ Print [ H, OnElementsOf Vol_FW, File StrCat[Dir,"map_H.pos" ]];
+
+ If(Flag_AnalysisType == 1) // current excitation
+ Print[ Rin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp0,"R=re(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempRin.dat"] ] ;
+ Print[ Xin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp1,"X=im(Z)"]{0}, Color "Ivory", File StrCat[Dir,"tempXin.dat"] ] ;
+ Print[ AbsZin, OnRegion Terminal, Format Table, SendToServer StrCat[pp2,"|Z|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsZin.dat"] ];
+ Print[ ArgZin, OnRegion Terminal, Format Table, SendToServer StrCat[pp3,"arg(Z) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgZin.dat"] ] ;
+ Else // voltage excitation
+ Print[ Gin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp0,"G=re(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempGin.dat"] ] ;
+ Print[ Bin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp1,"B=im(Y)"]{0}, Color "Ivory", File StrCat[Dir,"tempBin.dat"] ] ;
+ Print[ AbsYin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp2,"|Y|"]{0}, Color "Ivory", File StrCat[Dir,"tempAbsYin.dat"] ] ;
+ Print[ ArgYin, OnRegion Cut1H, Format Table, SendToServer StrCat[pp3,"arg(Y) [rad]"]{0}, Color "Ivory", File StrCat[Dir,"tempArgYin.dat"] ] ;
+ EndIf
+ Call Change_post_options;
+ }
+ }
+ EndIf
+}
diff --git a/ECE_full_wave_E-H_monopole/monopole_data.pro b/ECE_full_wave_E-H_monopole/monopole_data.pro
new file mode 100644
index 0000000000000000000000000000000000000000..ace76b6332bc782cd81ae643dc1728c70771836f
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/monopole_data.pro
@@ -0,0 +1,197 @@
+// Monopole
+// Geometries in separate .geo files
+//
+
+DefineConstant [ // allowing an external call (-set_number), e.g. from Matlab, python or a shell
+
+ // 1) spherical bnd, no PML, voltage excitation, ECE radiant, 2.5D, formulation in E, 1st order elements
+ BND_SHAPE = 2 // 1 for cylindrical, 2 for spherical, 2 for capsular
+ USE_PML = 0 // 1 for yes, 0 for no (radiant ECE)
+ TERMINAL_EXC = 0 // 0 for voltage excitation, 1 for current excitation
+ TYPE_BC_FictitiousBnd = 2 // 0 for non-terminal, 1 for ground, 2 for radiant (Silver-Muller)
+ FLAG_MODEL = 0 // 0 for 2.5 D, 1 for 3D
+ FORMULATION = 0 // 0 for E, 1 for H
+ FE_ORDER = 1 // 1 for order 1, 2 for order 2
+ //
+ FMIN = 0.5e9 // min freq [Hz]
+ FMAX = 4.5e9 // max freq [Hz]
+ NOP_FREQ = 50 // no of computed points
+ USE_FINE_MESH = 1 //1 - some very fine mesh coded in the geo file, 0 - only use the info below about nbLambda and nbDelta
+ NB_DELTA = 1 // no of elements per skin depth in the monopole, if USE_FINE_MESH = 1, otherwise not used
+ NOP_perLAMBDA_AIR = 5 // no of elements per minimum lambda in air; near the monopole it is finer, see the mesh description in the geo file
+ SCALE_CL = 1 // scale the characteristic length
+];
+
+// Geometry
+a = 1.52e-3; // monopole radius [m]
+b = 2.3*a; // minimum radius of the ground terminal
+L = 50e-3; // monopole length [m]
+rb = 4*L; // position of the fictitious boundary
+xb = rb; // needed for the cylindrical bnd and the capsular bnd
+yb = rb; // needed for the cylindrical bnd (not needed for the capsular bnd where yb = xb + L)
+zb = rb; // needed in 3D
+
+// ===============================
+// value of d indicates how the feeding is done:
+// d = 0 means no gap
+// d > 0 means a gap, the lenght of the gap is d
+// d < 0 is the coax case, -d is the length of the coax that feeds the monopole
+
+// uncomment the appropriate line below
+
+// no gap
+DefineConstant[
+ d = { 0, Choices{-4*a="-4*a", -2*a="-2*a", -a="-a", 0="0", a/2="a/2", a="a", 2*a="2*a"},
+ Name "{Geometrical Param./0Choose type of gap",
+ Help "if positive: gap with hard source feeding; if negative: coaxial with length -d",
+ Highlight "Red", Visible 1}
+
+ Flag_bndShape = {BND_SHAPE,
+ Choices{1="Cylindrical",2="Spherical", 3="Capsular"},
+ Name "{Geometrical Param./1Choose fictitious bnd shape",
+ Highlight "Tomato",Visible (!d)}
+
+ _use_pml = {USE_PML,
+ Choices{0,1}, Name "{Geometrical Param./2Use PML?",
+ Highlight "Pink",Visible (!d)}
+
+ Flag_AnalysisType = {TERMINAL_EXC,
+ Choices{0="voltage excitation", 1="current excitation"}, Highlight "ForestGreen",
+ Name "{Feed/Choose excitation type", ServerAction Str["Reset","GetDP/1ResolutionChoices"]}
+
+ Flag_AnalysisTypeExtBnd = {TYPE_BC_FictitiousBnd,
+ Choices{0="non-terminal",1="ground (PEC)", 2="Radiant ECE"},
+ Name "{Geometrical Param./3Choose BC backing the PML layer",
+ Highlight "Tomato",ServerAction Str["Reset", "Geometrical Param."],Visible (_use_pml)}
+
+ Flag_3Dmodel = {FLAG_MODEL,
+ Choices {0,1},
+ Name "{Geometrical Param./0Use 3D model?", Highlight "Pink",Visible (!d)}
+ ang_section_degrees = {360, Max 360, Name "{Geometrical Param./Wedge angle", Units "°", Visible Flag_3Dmodel,
+ Highlight "Pink"} // => cuts not be added for constraints
+
+ Flag_AnalysisTypeFormulation = {FORMULATION,
+ Choices{0="E-formulation", 1="H-formulation"}, Highlight "ForestGreen",
+ Name "{FE Param./Choose formulation", ServerAction Str["Reset","GetDP/1ResolutionChoices"]},
+
+ FE_Order = {FE_ORDER, Choices{1="1st order elements", 2="2nd order elements (hierarchical)"}, Highlight "ForestGreen",
+ Name "{FE Param./Choose order", ServerAction Str["Reset","GetDP/1ResolutionChoices"]}
+];
+
+// equivalent to previous If-Else-EndIf
+_use_radiantECE = (Flag_AnalysisTypeExtBnd >= 2);
+Flag_SilverMuller = (Flag_AnalysisTypeExtBnd == 2);
+Flag_PEC = (Flag_AnalysisTypeExtBnd == 3);
+
+modelDim = (!Flag_3Dmodel)? 2: 3;
+Flag_Axi = (!Flag_3Dmodel)? 1: 0; // If flag set to 0, problem is either 2D planar or 3D
+h2Ddepth = (!Flag_3Dmodel)? 2*Pi:360/ang_section_degrees;
+
+// Frequencies
+fmin = FMIN;
+fmax = FMAX;
+nop = NOP_FREQ;
+//freqs() = LogSpace[Log10[fmin],Log10[fmax],nop];
+freqs() = LinSpace[fmin,fmax,nop];
+
+cGUI2= "{Monopole antenna1";
+mGeoF = StrCat[cGUI2,"Frequency information/0"];
+colorMGeoF = "Orange";
+close_menu = 1;
+DefineConstant[
+ _use_freq_loop = {0, Choices{0,1}, Name StrCat[mGeoF, "0Loop on frequencies?"],
+ ServerAction Str["Reset", StrCat[mGeoF, "0Working Frequency"]], Highlight Str[colorMGeoF]}
+ Freq = {freqs(0), Choices{freqs()},
+ Loop _use_freq_loop, Name StrCat[mGeoF, "0Working Frequency"],
+ Units "Hz", Highlight Str[colorMGeoF], Closed !close_menu }
+
+ PmlDelta = {L/2, Visible _use_pml, // Always defined, but used only if PML, possible freq dependance
+ Name StrCat[mGeoF, "0PML thickness"], Units "m" }
+
+ _use_fineMesh = {USE_FINE_MESH, Choices{0="coarse", 1="fine"}, Highlight "ForestGreen",
+ Name "{FE Mesh/Choose fineness", ServerAction Str["Reset","GetDP/1ResolutionChoices"]}
+];
+
+nopPerLambdaAir = NOP_perLAMBDA_AIR;
+s = SCALE_CL; // a factor that multiplies each characteristic length to easily refine the mesh
+nbDelta = (!_use_fineMesh) ? NB_DELTA :1;
+
+// ================================= Material
+// Monopole
+epsrMonopole = 1;
+murMonopole = 1;
+sigmaMonopole = 6.6e7;
+// Air
+epsrAir = 1;
+murAir = 1;
+sigmaAir = 0;
+// Coax conductor - used only if d < 0
+epsrCoaxCond = 1;
+murCoaxCond = 1;
+sigmaCoaxCond = 6.6e7;
+// Coax dielectric - used only if d < 0
+epsrCoaxDiel = 1;
+murCoaxDiel = 1;
+sigmaCoaxDiel = 1e-3;
+
+eps0 = 8.854187818e-12; // [F/m] - absolute permitivity of vacuum
+mu0 = 4.e-7 * Pi; // [H/m] - absolute permeability of vacuum
+nu0 = 1/mu0;
+
+// epsilon and nu for all the materials
+epsMonopole = eps0*epsrMonopole;
+nuMonopole = 1/(mu0*murMonopole);
+
+epsAir = eps0*epsrAir;
+nuAir = 1/(mu0*murAir);
+
+epsCoaxCond = eps0*epsrCoaxCond;
+nuCoaxCond = 1/(mu0*murCoaxCond);
+
+epsCoaxDiel = eps0*epsrCoaxDiel;
+nuCoaxDiel = 1/(mu0*murCoaxDiel);
+
+c0 = 3e8; // [m/s] speed of light in vacuum
+lambda = c0/Freq;
+k0 = 2*Pi/lambda; // Wave number
+
+
+// ================================= Some quantities output to GUI
+cGUI8= "{Monopole antenna/0";
+If (Flag_AnalysisTypeFormulation == 0) // E
+ pp0 = StrCat[cGUI8,"Output-e/0"];
+ pp1 = StrCat[cGUI8,"Output-e/1"];
+ pp2 = StrCat[cGUI8,"Output-e/2"];
+ pp3 = StrCat[cGUI8,"Output-e/3"];
+Else
+ pp0 = StrCat[cGUI8,"Output-h/0"];
+ pp1 = StrCat[cGUI8,"Output-h/1"];
+ pp2 = StrCat[cGUI8,"Output-h/2"];
+ pp3 = StrCat[cGUI8,"Output-h/3"];
+EndIf
+
+//==================================
+// Indexes of physical entities
+// Surfaces for 2.5D and volumes for 3D
+MONOPOLE = 1000;
+AIR = 3000;
+PML = 3100 ;
+
+// boundaries - lines for 2.5D and surfaces for 3D
+TERMINAL = 2000 ;
+GROUND = 2100 ;
+CUT = 2300 ;
+
+SURFAIRINF = 3333 ;
+AXIS = 4444 ;
+
+// used in 3D
+NONTERM = 4000;
+CUT1H = NONTERM+1;
+
+//=======================================================================
+// Indexes of physical entities when Gap is non zero or Coax feeding is used
+// Surfaces
+COAXCONDUCTOR = 1001; // for d < 0
+COAXDIELECTRIC = 3001;
+
diff --git a/ECE_full_wave_E-H_monopole/pro_problemIndependent/FullWave_E_eceRadiant_SISO_secondOrder.pro b/ECE_full_wave_E-H_monopole/pro_problemIndependent/FullWave_E_eceRadiant_SISO_secondOrder.pro
new file mode 100644
index 0000000000000000000000000000000000000000..da579e0c6e233285c7de515287a8834ec03f7520
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/pro_problemIndependent/FullWave_E_eceRadiant_SISO_secondOrder.pro
@@ -0,0 +1,112 @@
+/* FullWave_E_ece_SISO.pro
+
+ Problem independent pro file
+ for passive, linear domains in Full Wave with radiant ECE boundary conditions.
+
+ It uses a formulation with vec{E} strictly inside the domain and V strictly on the boundary.
+ The domain is simply connected, it has only electric terminals.
+ The terminals can be voltage or current excited.
+ The material inside is linear from all points of view: electric, magnetic, conduction
+ It is compulsory that there exists a ground terminal (its voltage is set to zero).
+ The post-operation assumes that there is one excited (active) terminal (excitation equal to 1 - it can be in
+ voltage or current excited terminal).
+ Frequency domain simulations.
+*/
+
+Include "Jacobian_and_Integration.pro"
+Include "only_FunctionSpace_and_Formulation_FullWave_E_eceRadiant_secondOrder.pro"
+
+Resolution {
+ { Name FullWave_E_ece;
+ System {
+ { Name Sys_Ele; NameOfFormulation FullWave_E_ece; Type Complex; Frequency Freq;}
+ }
+ Operation {
+ CreateDir[StrCat[modelDir,Dir]];
+ SetTime[Freq]; // usefull to save (and print) the current frequency
+ //GenerateSeparate[Sys_Ele]; // separate matrices
+ Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
+ //Print[Sys_Ele]; // save in matlab format or change the format when running getdp
+ }
+ }
+}
+
+PostProcessing {
+
+ { Name FW_E_ece; NameOfFormulation FullWave_E_ece;
+ Quantity {
+ { Name E; Value {Local { [ {e} ]; In Vol_FW; Jacobian Vol; }}} // complex
+ { Name rmsE; Value { Local { [ Norm[{e}] ]; In Vol_FW; Jacobian Vol; }}} // complex
+ { Name J; Value { Local { [ sigma[]*{e} ]; In Vol_FW; Jacobian Vol; }}} // complex
+ { Name rmsJ; Value { Local { [ Norm[sigma[]*{e}] ]; In Vol_FW; Jacobian Vol; }}} // real
+ { Name gradV; Value { Local { [ {dv} ]; In Vol_FW; Jacobian Vol; }}} // complex
+
+ { Name Vterminals;
+ Value {
+ // This is a global variable, you do not need to indicate the Jacobian
+ // It is only defined at the terminal and it is only one value for the whole line (surface in 3D)
+ // If the support you give here is wider, you will get 0 where it is not defined
+ Term { [ -{V} ]; In Sur_Terminals_FWece; } // => unknowns only in Sur_Terminals_FWece, outside it will be zero
+ }
+ }
+
+ { Name Rin; Value { Term { [ Re[-{V}] ]; In Sur_Terminals_FWece; }}}
+ { Name Xin; Value { Term { [ Im[-{V}] ]; In Sur_Terminals_FWece; }}}
+ { Name AbsZin; Value { Term { [Norm[-{V}] ]; In Sur_Terminals_FWece; }}}
+ { Name ArgZin; Value { Term { [Atan2[ Im[-{V}], Re[-{V}]] ]; In Sur_Terminals_FWece; }}}
+
+ { Name I;
+ Value { // h2Ddepth is the depth for 2D problems; =1 for 3D problems
+ Term { [ -{I}*h2Ddepth ]; In Sur_Terminals_FWece; }
+ }
+ }
+
+ // Yin = I/V = Gin + j Bin
+ { Name Yin; Value { Term { [ -{I}*h2Ddepth * Conj[-{V}]/SquNorm[{V}] ]; In Sur_Terminals_FWece; }}}
+ { Name AbsYin; Value { Term { [ Norm[-{I}*h2Ddepth * Conj[-{V}]/SquNorm[{V}]] ]; In Sur_Terminals_FWece; }}}
+ { Name ArgYin; Value { Term { [ Atan2[ Im[-{I}*h2Ddepth * Conj[-{V}]/SquNorm[{V}]], Re[ -{I}*h2Ddepth * Conj[-{V}]/SquNorm[{V}]] ] ]; In Sur_Terminals_FWece; }}}
+
+ { Name Gin; Value { Term { [Re[ -{I}*h2Ddepth ]]; In Sur_Terminals_FWece; }}}
+ { Name Bin; Value { Term { [Im[ -{I}*h2Ddepth ]]; In Sur_Terminals_FWece; }}}
+ { Name AbsYin_; Value { Term { [ Norm[-{I}*h2Ddepth ] ]; In Sur_Terminals_FWece; }}}
+ { Name ArgYin_; Value { Term { [ Atan2[ Im[-{I}*h2Ddepth ], Re[-{I}*h2Ddepth ]] ]; In Sur_Terminals_FWece; }}}
+
+ { Name VnotTerminals;
+ Value {
+ Local { [ -{dInv dv} ]; In Vol_FW; Jacobian Vol; }
+ Local { [ -{dInv dvf} ]; In Vol_FW; Jacobian Vol; }
+ }
+ }
+
+ { Name B;
+ Value {
+ Local { [ {d e}/(2*Pi*Freq*Complex[0,1]) ]; In Vol_FW; Jacobian Vol; }
+ }
+ }
+
+ { Name H;
+ Value {
+ Local { [ -nu[]*{d e}/(2*Pi*Freq*Complex[0,1]) ]; In Vol_FW; Jacobian Vol; }
+ }
+ }
+
+ { Name rmsH;
+ Value {
+ Local { [Norm[ -nu[]*{d e}/(2*Pi*Freq*Complex[0,1]) ]]; In Vol_FW; Jacobian Vol; }
+ }
+ }
+
+ { Name nVec;
+ Value {
+ Local { [ XYZ[]/Norm[XYZ[]] ]; In Vol_FW; Jacobian Vol; }
+ }
+ }
+
+ { Name tVec;
+ Value {
+ Local { [ (XYZ[]/Norm[XYZ[]]) /\ Vector[0,0,1] ]; In Vol_FW; Jacobian Vol; }
+ }
+ }
+ }
+ }
+ }
diff --git a/ECE_full_wave_E-H_monopole/pro_problemIndependent/FullWave_H_eceRadiant_SISO_cplx_2D_and_AXI_secondOrder.pro b/ECE_full_wave_E-H_monopole/pro_problemIndependent/FullWave_H_eceRadiant_SISO_cplx_2D_and_AXI_secondOrder.pro
new file mode 100644
index 0000000000000000000000000000000000000000..d6640ff3d5a8a7c2afc05de68078adb8ce681044
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/pro_problemIndependent/FullWave_H_eceRadiant_SISO_cplx_2D_and_AXI_secondOrder.pro
@@ -0,0 +1,113 @@
+/* FullWave_H_ece_SISO_cplx_2D.pro/
+
+ Problem independent pro file
+ for passive, linear domains in Full Wave with radiant ECE boundary conditions.
+
+ It uses a formulation with vec{H} strictly inside the domain and V strictly on the boundary.
+ The domain is simply connected, it has only electric terminals.
+ The terminals can be voltage or current excited.
+ The material inside is linear from all points of view: electric, magnetic, conduction
+ It is compulsory that there exists a ground terminal (its voltage is set to zero).
+ The post-operation assumes that there is one excited (active) terminal (excitation equal to 1 - it can be in
+ voltage or current excited terminal).
+ Frequency domain simulations.
+*/
+
+Include "Jacobian_and_Integration.pro"
+Include "only_FunctionSpace_and_Formulation_FullWave_H_eceRadiant_cplx_2D_and_AXI_secondOrder.pro"
+
+Resolution {
+ { Name FullWave_H_ece;
+ System {
+ { Name Sys_Ele; NameOfFormulation FullWave_H_ece; Type Complex; Frequency Freq;}
+ }
+ Operation {
+ CreateDir[StrCat[modelDir,Dir]];
+ SetTime[Freq]; // usefull to save (and print) the current frequency
+ //GenerateSeparate[Sys_Ele]; // separate matrices
+ Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
+ //Print[Sys_Ele]; // save in matlab format or change the format when running getdp
+ }
+ }
+}
+
+
+PostProcessing {
+
+ { Name FW_H_ece; NameOfFormulation FullWave_H_ece;
+ Quantity {
+ { Name H; Value { Local { [ {h} ]; In Vol_FW; Jacobian Vol; }}}
+ { Name Hz; Value { Local { [ CompZ[{h}] ]; In Vol_FW; Jacobian Vol; }}}
+ { Name rmsH; Value { Local { [ Norm[{h}] ]; In Vol_FW; Jacobian Vol; }}}
+
+
+ { Name E; Value { Local { [ {d h}/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[]) ]; In Vol_FW; Jacobian Vol; }}}
+ { Name Ey; Value { Local { [ CompY[{d h}/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[])] ]; In Vol_FW; Jacobian Vol; }}}
+
+ { Name J; Value { Local { [ {d h}*sigma[]/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[]) ]; In Vol_FW; Jacobian Vol; }}}
+ { Name rmsJ; Value { Local { [ Norm[{d h}*sigma[]/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[])] ]; In Vol_FW; Jacobian Vol; }}}
+
+
+ { Name Vterminals; Value { Term { [ -{V} ]; In Cut1H; }}}
+
+ { Name Rin; Value { Term { [ Re[-{V}] ]; In Cut1H; }}}
+ { Name Xin; Value { Term { [ Im[-{V}] ]; In Cut1H; }}}
+ { Name AbsZin; Value { Term { [ Norm[ -{V} ] ]; In Cut1H; }}}
+ { Name ArgZin; Value { Term { [Atan2[Im[-{V} ],Re[ -{V} ]]]; In Cut1H; }}}
+
+ { Name I;
+ Value { // h is the depth for 2D problems, and it has to be 1 for 3D problems
+ If(Flag_Axi == 0)
+ Term { [ -2*{I}*h2Ddepth ]; In Cut1H; }
+ Else // 2D axisymmetric
+ Term { [ -{I}*h2Ddepth ]; In Cut1H; }
+ EndIf
+ }
+ }
+
+ { Name Gin;
+ Value {
+ If(Flag_Axi == 0) // 2D
+ Term { [Re[ -2*{I}*h2Ddepth ]]; In Cut1H; }
+ Else
+ Term { [Re[ -{I}*h2Ddepth ]]; In Cut1H; }
+ EndIf
+ }
+ }
+
+ { Name Bin;
+ Value {
+ If(Flag_Axi == 0) // 2D
+ Term { [Im[ -2*{I}*h2Ddepth ]]; In Cut1H; }
+ Else
+ Term { [Im[ -{I}*h2Ddepth ]]; In Cut1H; }
+ EndIf
+ }
+ }
+
+ { Name AbsYin;
+ Value {
+ If(Flag_Axi == 0) // 2D
+ Term { [ Norm[ -2*{I}*h2Ddepth ] ]; In Cut1H; }
+ Else
+ Term { [ Norm[ -{I}*h2Ddepth ] ]; In Cut1H; }
+ EndIf
+ }
+ }
+
+ { Name ArgYin;
+ Value {
+ If(Flag_Axi == 0) // 2D
+ Term { [Atan2[Im[ -2*{I}*h2Ddepth ],Re[ -2*{I}*h2Ddepth ]]]; In Cut1H; }
+ Else
+ Term { [Atan2[Im[ -{I}*h2Ddepth ],Re[ -{I}*h2Ddepth ]]]; In Cut1H; }
+ EndIf
+ }
+ }
+ }
+ }
+ }
+
+
+
+
diff --git a/ECE_full_wave_E-H_monopole/pro_problemIndependent/FullWave_H_eceRadiant_SISO_cplx_3D_secondOrder.pro b/ECE_full_wave_E-H_monopole/pro_problemIndependent/FullWave_H_eceRadiant_SISO_cplx_3D_secondOrder.pro
new file mode 100644
index 0000000000000000000000000000000000000000..5a2154bcfa5a23e9cb9188850e7d604edc98bf7e
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/pro_problemIndependent/FullWave_H_eceRadiant_SISO_cplx_3D_secondOrder.pro
@@ -0,0 +1,65 @@
+/* FullWave_H_eceRadiant_SISO_cplx_3D_secondOrder.pro
+
+ Problem independent pro file
+ for passive, linear domains in Full Wave with ECE boundary conditions.
+
+ It uses a formulation with vec{H} strictly inside the domain and V strictly on the boundary.
+ The domain is simply connected, it has only electric terminals.
+ The terminals can be voltage or current excited.
+ The material inside is linear from all points of view: electric, magnetic, conduction
+ It is compulsory that there exists a ground terminal (its voltage is set to zero).
+ The post-operation assumes that there is one excited (active) terminal (excitation equal to 1 - it can be in
+ voltage or current excited terminal).
+ Frequency domain simulations.
+*/
+
+Include "Jacobian_and_Integration.pro"
+Include "only_FunctionSpace_and_Formulation_FullWave_H_eceRadiant_cplx_3D_secondOrder.pro"
+
+Resolution {
+ { Name FullWave_H_ece;
+ System {
+ { Name Sys_Ele; NameOfFormulation FullWave_H_ece; Type Complex; Frequency Freq;}
+ }
+ Operation {
+ CreateDir[StrCat[modelDir,Dir]];
+ SetTime[Freq]; // usefull to save (and print) the current frequency
+ //GenerateSeparate[Sys_Ele]; // separate matrices
+ Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
+ //Print[Sys_Ele]; // save in matlab format or change the format when running getdp
+ }
+ }
+}
+
+
+PostProcessing {
+
+ { Name FW_H_ece; NameOfFormulation FullWave_H_ece;
+ Quantity {
+ { Name H; Value { Local { [ {h} ]; In Vol_FW; Jacobian Vol; }}} // complex
+ { Name Hz ; Value { Local{ [ CompZ[{h}] ]; In Vol_FW ; Jacobian Vol ;}}}
+ { Name rmsH; Value { Local { [ Norm[{h}] ]; In Vol_FW; Jacobian Vol; }}}
+
+
+ { Name E; Value { Local { [ {d h}/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[]) ]; In Vol_FW; Jacobian Vol; }}}
+ { Name Ey; Value { Local { [ CompY[{d h}/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[])] ]; In Vol_FW; Jacobian Vol; }}}
+ { Name J; Value { Local { [ {d h}*sigma[]/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[]) ]; In Vol_FW; Jacobian Vol; }}}
+ { Name rmsJ; Value { Local { [ Norm[{d h}*sigma[]/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[])] ]; In Vol_FW; Jacobian Vol; }}}
+
+
+ { Name Vterminals; Value { Term { [ -{V} ]; In Cut1H; }}}
+
+ { Name Rin; Value {Term { [ Re[-{V}] ]; In Cut1H; }}}
+ { Name Xin; Value {Term { [ Im[-{V}] ]; In Cut1H; }}}
+ { Name AbsZin; Value { Term { [ Norm[ -{V} ] ]; In Cut1H; }}}
+ { Name ArgZin; Value { Term { [ Atan2[Im[ -{V} ], Re[ -{V} ]]]; In Cut1H; }}}
+
+ { Name I; Value { Term { [ -{I}*h2Ddepth ]; In Cut1H; }}} // h2depth = 1 for 3D problems
+
+ { Name Gin; Value { Term { [Re[ -{I}*h2Ddepth ]]; In Cut1H; }}}
+ { Name Bin; Value { Term { [Im[ -{I}*h2Ddepth ]]; In Cut1H; }}}
+ { Name AbsYin; Value { Term { [ Norm[ -{I}*h2Ddepth ]]; In Cut1H; }}}
+ { Name ArgYin; Value { Term { [Atan2[Im[ -{I}*h2Ddepth ], Re[ -{I}*h2Ddepth ]]]; In Cut1H; }}}
+ }
+ }
+ }
diff --git a/ECE_full_wave_E-H_monopole/pro_problemIndependent/Jacobian_and_Integration.pro b/ECE_full_wave_E-H_monopole/pro_problemIndependent/Jacobian_and_Integration.pro
new file mode 100644
index 0000000000000000000000000000000000000000..e67f5f8a6e6c5a1b6c12954042f2e6d13e4fb5fc
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/pro_problemIndependent/Jacobian_and_Integration.pro
@@ -0,0 +1,68 @@
+Jacobian {
+ { Name Vol;
+ Case {
+ If(Flag_Axi) // 2D axisymetric problems
+ { Region All; Jacobian VolAxiSqu; }
+ Else
+ { Region All; Jacobian Vol; }
+ EndIf
+ }
+ }
+ { Name Sur;
+ Case {
+ If(Flag_Axi) // 2D axisymetric problems
+ { Region All; Jacobian SurAxi; }
+ Else
+ { Region All; Jacobian Sur; }
+ EndIf
+ }
+ }
+}
+
+Integration {
+ { Name Int;
+ If(FE_Order !=2) // Adapt gauss integration points with higher order
+ Case {
+ { Type Gauss;
+ Case {
+ { GeoElement Point; NumberOfPoints 1; }
+ { GeoElement Line; NumberOfPoints 3; }
+ { GeoElement Triangle; NumberOfPoints 3; }
+ { GeoElement Quadrangle; NumberOfPoints 4; }
+ { GeoElement Tetrahedron; NumberOfPoints 4; }
+ { GeoElement Hexahedron; NumberOfPoints 6; }
+ { GeoElement Prism; NumberOfPoints 9; }
+ { GeoElement Pyramid; NumberOfPoints 8; }
+
+ // 2nd order with geometry
+ { GeoElement Line2; NumberOfPoints 4; }
+ { GeoElement Triangle2; NumberOfPoints 12; }
+ { GeoElement Quadrangle2; NumberOfPoints 12; }
+ { GeoElement Tetrahedron2; NumberOfPoints 17; }
+ }
+ }
+ }
+ Else
+ Case {
+ { Type Gauss;
+ Case {
+ { GeoElement Point; NumberOfPoints 1; }
+ { GeoElement Line; NumberOfPoints 10; }
+ { GeoElement Triangle; NumberOfPoints 16 ; }
+ { GeoElement Quadrangle; NumberOfPoints 7; }
+ { GeoElement Tetrahedron; NumberOfPoints 15; }
+ { GeoElement Hexahedron; NumberOfPoints 6; }
+ { GeoElement Prism; NumberOfPoints 9; }
+ { GeoElement Pyramid; NumberOfPoints 8; }
+
+ // 2nd order with geometry
+ { GeoElement Line2; NumberOfPoints 4; }
+ { GeoElement Triangle2; NumberOfPoints 12; }
+ { GeoElement Quadrangle2; NumberOfPoints 12; }
+ { GeoElement Tetrahedron2; NumberOfPoints 17; }
+ }
+ }
+ }
+ EndIf
+ }
+}
diff --git a/ECE_full_wave_E-H_monopole/pro_problemIndependent/only_FunctionSpace_and_Formulation_FullWave_E_eceRadiant_secondOrder.pro b/ECE_full_wave_E-H_monopole/pro_problemIndependent/only_FunctionSpace_and_Formulation_FullWave_E_eceRadiant_secondOrder.pro
new file mode 100644
index 0000000000000000000000000000000000000000..4a5bdb86f971c7335e5176112ffeef3fbc6496d3
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/pro_problemIndependent/only_FunctionSpace_and_Formulation_FullWave_E_eceRadiant_secondOrder.pro
@@ -0,0 +1,112 @@
+FunctionSpace {
+
+ { Name Hcurl_E; Type Form1;
+ BasisFunction {
+ { Name se; NameOfCoef ee; Function BF_Edge;
+ Support Dom_FW ; Entity EdgesOf[All, Not Sur_FW]; }
+ // additional basis functions for 2nd order interpolation (hierachical)
+ If(FE_Order == 2)
+ { Name se2; NameOfCoef ee2; Function BF_Edge_2E;
+ Support Dom_FW; Entity EdgesOf[All, Not Sur_FW]; }
+ // Reduced 2nd order
+ { Name se3a ; NameOfCoef ee3a ; Function BF_Edge_3F_a ;
+ Support Dom_FW ; Entity FacetsOf[All, Not Sur_FW] ; }
+ { Name se3b ; NameOfCoef ee3b ; Function BF_Edge_3F_b ;
+ Support Dom_FW ; Entity FacetsOf[ All, Not Sur_FW] ; }
+ // 2nd order
+ { Name se4; NameOfCoef ee4; Function BF_Edge_4E;
+ Support Dom_FW ; Entity EdgesOf[All, Not Sur_FW]; }
+ EndIf
+ If (_use_radiantECE)
+ { Name sn; NameOfCoef vn; Function BF_GradNode;
+ Support Dom_FW ; Entity NodesOf[Sur_FW, Not Sur_Terminals_FWeceAndInf]; }
+ If(FE_Order == 2)
+ // GradNode BFs
+ { Name sn2; NameOfCoef vn2; Function BF_GradNode_2E;
+ Support Dom_FW ; Entity EdgesOf[Sur_FW, Not Sur_Terminals_FWeceAndInf]; }
+ EndIf
+ Else
+ { Name sn; NameOfCoef vn; Function BF_GradNode;
+ Support Dom_FW ; Entity NodesOf[Sur_FW, Not Sur_Terminals_FWece]; }
+ If(FE_Order == 2)
+ // GradNode BFs
+ { Name sn2; NameOfCoef vn2; Function BF_GradNode_2E;
+ Support Dom_FW ; Entity EdgesOf[Sur_FW, Not Sur_Terminals_FWece]; }
+ EndIf
+ EndIf
+
+ { Name sf; NameOfCoef vf; Function BF_GradGroupOfNodes;
+ Support Dom_FW ; Entity GroupsOfNodesOf[Sur_Terminals_FWece]; }
+ If (_use_radiantECE)
+ { Name seInf; NameOfCoef eeInf; Function BF_Edge;
+ Support Dom_FW ; Entity EdgesOf[SurfAirInf]; }
+ If (1==0) // (FE_Order == 2)
+ { Name se2Inf; NameOfCoef ee2Inf; Function BF_Edge_2E;
+ Support Dom_FW; Entity EdgesOf[SurfAirInf]; }
+ // Reduced 2nd order
+ { Name se3aInf ; NameOfCoef ee3aInf ; Function BF_Edge_3F_a ;
+ Support Dom_FW ; Entity FacetsOf[SurfAirInf] ; }
+ { Name se3bInf ; NameOfCoef ee3bInf ; Function BF_Edge_3F_b ;
+ Support Dom_FW ; Entity FacetsOf[SurfAirInf] ; }
+ // 2nd order
+ { Name se4Inf; NameOfCoef ee4Inf; Function BF_Edge_4E;
+ Support Dom_FW ; Entity EdgesOf[SurfAirInf]; }
+ EndIf
+ EndIf
+ }
+ GlobalQuantity {
+ { Name TerminalPotential; Type AliasOf ; NameOfCoef vf; }
+ { Name TerminalCurrent ; Type AssociatedWith; NameOfCoef vf; }
+ }
+
+ SubSpace {
+ { Name dv ; NameOfBasisFunction {sn}; } // Subspace, it maybe use in equations or post-pro
+ { Name dvf ; NameOfBasisFunction {sf}; }
+ }
+
+ Constraint {
+ { NameOfCoef TerminalPotential; EntityType GroupsOfNodesOf;
+ NameOfConstraint SetTerminalPotential; }
+ { NameOfCoef TerminalCurrent; EntityType GroupsOfNodesOf;
+ NameOfConstraint SetTerminalCurrent; }
+ If (_use_radiantECE)
+ If (Flag_SilverMuller==0) // classic, PEC
+ If (Flag_PEC == 1)
+ { NameOfCoef eeInf; EntityType EdgesOf; //ee are voltages along edges
+ NameOfConstraint ElectricVoltages;} // impose voltages along edges with known Et
+ EndIf
+ EndIf
+ EndIf
+ }
+ }
+}
+
+Formulation {
+
+ { Name FullWave_E_ece; Type FemEquation;
+ Quantity {
+ { Name e; Type Local; NameOfSpace Hcurl_E; }
+ { Name dv; Type Local; NameOfSpace Hcurl_E[dv]; } // Just for post-processing issues
+ { Name dvf; Type Local; NameOfSpace Hcurl_E[dvf]; }
+ { Name V; Type Global; NameOfSpace Hcurl_E[TerminalPotential]; }
+ { Name I; Type Global; NameOfSpace Hcurl_E[TerminalCurrent]; }
+ }
+ Equation {
+ // \int_D curl(\vec{E}) \cdot curl(\vec{e}) dv
+ Galerkin { [ nu[] * Dof{d e} , {d e} ]; In Vol_FW; Jacobian Vol; Integration Int; }
+
+ // \int_D j*\omega*(\sigma + j*\omega*\epsilon) \vec{E} \cdot \vec{e} dv
+ Galerkin { DtDof [ sigma[] * Dof{e} , {e} ]; In Vol_FW; Jacobian Vol; Integration Int; }
+ Galerkin { DtDtDof [ epsilon[] * Dof{e} , {e} ]; In Vol_FW; Jacobian Vol; Integration Int; }
+
+ // j*\omega*sum (vk Ik); for k - all terminals so that the currents through the terminals will be computed as well
+ GlobalTerm {DtDof [ -Dof{I} , {V} ]; In Sur_Terminals_FWece; }
+
+ If (_use_radiantECE && Flag_SilverMuller)
+ Galerkin { DtDof [ Sqrt[epsilonNu[]] * ( Normal[] /\ Dof{e} ) /\ Normal[] , {e} ];
+ In SurfAirInf; Jacobian Sur; Integration Int; }
+ EndIf
+
+ }
+ }
+}
diff --git a/ECE_full_wave_E-H_monopole/pro_problemIndependent/only_FunctionSpace_and_Formulation_FullWave_H_eceRadiant_cplx_2D_and_AXI_secondOrder.pro b/ECE_full_wave_E-H_monopole/pro_problemIndependent/only_FunctionSpace_and_Formulation_FullWave_H_eceRadiant_cplx_2D_and_AXI_secondOrder.pro
new file mode 100644
index 0000000000000000000000000000000000000000..0ea421f5bfe8dd10726b7f0d012c3f04bcee38bc
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/pro_problemIndependent/only_FunctionSpace_and_Formulation_FullWave_H_eceRadiant_cplx_2D_and_AXI_secondOrder.pro
@@ -0,0 +1,98 @@
+//=========================================================
+// Function spaces
+//=========================================================
+
+FunctionSpace {
+
+ { Name Hcurl_H; Type Form1P; // Hspace for 2D problems with H perpendicular to the plane
+ BasisFunction {
+
+ If (_use_radiantECE)
+ { Name se; NameOfCoef he; Function BF_PerpendicularEdge;
+ Support Dom_FW; Entity NodesOf[All, Not Cut1H]; }
+ { Name sc; NameOfCoef ic; Function BF_GroupOfPerpendicularEdges;
+ Support Dom_FW; Entity GroupsOfNodesOf[Cut1H]; }
+
+ If(FE_Order == 2) // hierarchical basis functions for 2nd order interpolation
+ { Name se2; NameOfCoef he2; Function BF_PerpendicularEdge_2E;
+ Support Dom_FW; Entity EdgesOf[All, Not Cut1H]; }
+ { Name sf2; NameOfCoef hf2; Function BF_PerpendicularEdge_2F;
+ Support Dom_FW; Entity FacetsOf[All, Not Cut1H]; }
+ EndIf
+ Else
+ { Name se; NameOfCoef he; Function BF_PerpendicularEdge;
+ Support Dom_FW; Entity NodesOf[All, Not BoundaryNotTerminal]; }
+ { Name sc; NameOfCoef ic; Function BF_GroupOfPerpendicularEdges;
+ Support Dom_FW; Entity GroupsOfNodesOf[BoundaryNotTerminal]; }
+
+ If(FE_Order == 2) // hierarchical basis functions for 2nd order interpolation
+ { Name se2; NameOfCoef he2; Function BF_PerpendicularEdge_2E;
+ Support Dom_FW; Entity EdgesOf[All, Not BoundaryNotTerminal]; }
+ { Name sf2; NameOfCoef hf2; Function BF_PerpendicularEdge_2F;
+ Support Dom_FW; Entity FacetsOf[All, Not BoundaryNotTerminal]; }
+ EndIf
+ EndIf
+ }
+
+ GlobalQuantity {
+ { Name TerminalCurrent; Type AliasOf ; NameOfCoef ic; }
+ { Name TerminalPotential; Type AssociatedWith; NameOfCoef ic; }
+ }
+
+ Constraint {
+ { NameOfCoef TerminalPotential; EntityType Auto; NameOfConstraint SetTerminalPotentialH; }
+ { NameOfCoef TerminalCurrent; EntityType Auto; NameOfConstraint SetTerminalCurrentH; }
+
+ { NameOfCoef he; EntityType Auto; NameOfConstraint ImposeHonAxis; }
+ If(FE_Order==2)
+ { NameOfCoef he2; EntityType Auto; NameOfConstraint ImposeHonAxis; }
+ { NameOfCoef hf2; EntityType Auto; NameOfConstraint ImposeHonAxis; }
+ EndIf
+ If (_use_radiantECE && Flag_SilverMuller==0 && Flag_PEC == 0)
+ Printf(" ======================================================================== mg voltages set to zero on SurfInf");
+ { NameOfCoef he; EntityType EdgesOf; // he are magnetic voltages along edges
+ NameOfConstraint MagneticVoltages;} // impose voltages along edges with known Ht
+ EndIf
+ }
+ }
+}
+
+Function{
+ Omega = 2*Pi*Freq;
+ jOmega[] = Complex[0,1]*Omega;
+ jOmega2[] = -Omega^2;
+
+ tensSuma[] = sigma[]/jOmega[] + epsilon[];
+ alpha[] = 1/tensSuma[];
+ alphaIm[] = Complex[0,1]*Im[alpha[]];
+ alphaRe[] = Re[alpha[]];
+
+}
+
+Formulation {
+
+ { Name FullWave_H_ece; Type FemEquation;
+ Quantity {
+ { Name h; Type Local; NameOfSpace Hcurl_H; }
+ { Name V; Type Global; NameOfSpace Hcurl_H[TerminalPotential]; }
+ { Name I; Type Global; NameOfSpace Hcurl_H[TerminalCurrent]; }
+ }
+ Equation {
+ // multiplying with j omega
+ // \int_D \alpha * curl(\vec{H}) \cdot curl(\vec{h}) dv
+
+ Galerkin { [ alpha[] * Dof{d h} , {d h} ]; In Vol_FW; Jacobian Vol; Integration Int; }
+
+ // \int_D j*\omega*(j*\omega*\mu) \vec{H} \cdot \vec{h} dv
+ Galerkin { [ jOmega2[] * mu[] * Dof{h} , {h} ]; In Vol_FW; Jacobian Vol; Integration Int; }
+
+ // sum (Vk ik); for k - all terminals so that the terminal voltages will be computed as well
+ GlobalTerm { [ jOmega[]*Dof{V}, {I} ]; In Cut1H; }
+
+ If (_use_radiantECE && Flag_SilverMuller)
+ Galerkin { DtDof [ Sqrt[1/epsilonNu[]] * ( Normal[] /\ Dof{h} ) /\ Normal[] , {h} ];
+ In SurfAirInf; Jacobian Sur; Integration Int; }
+ EndIf
+ }
+ }
+}
diff --git a/ECE_full_wave_E-H_monopole/pro_problemIndependent/only_FunctionSpace_and_Formulation_FullWave_H_eceRadiant_cplx_3D_secondOrder.pro b/ECE_full_wave_E-H_monopole/pro_problemIndependent/only_FunctionSpace_and_Formulation_FullWave_H_eceRadiant_cplx_3D_secondOrder.pro
new file mode 100644
index 0000000000000000000000000000000000000000..1a8643391cfb05810c7d300e8943dd8404d33ba6
--- /dev/null
+++ b/ECE_full_wave_E-H_monopole/pro_problemIndependent/only_FunctionSpace_and_Formulation_FullWave_H_eceRadiant_cplx_3D_secondOrder.pro
@@ -0,0 +1,145 @@
+//=========================================================
+// Function spaces
+//=========================================================
+
+FunctionSpace {
+
+ { Name Hcurl_H; Type Form1; // Hspace for 3D
+ BasisFunction {
+
+ If (_use_radiantECE)
+ { Name se; NameOfCoef he; Function BF_Edge;
+ Support Dom_FW; Entity EdgesOf[All, Not NonTerm]; }
+
+ { Name sn; NameOfCoef phin; Function BF_GradNode;
+ Support Dom_FW;
+ Entity NodesOf[BoundaryNotTerminal, Not SurfAirInf]; }
+
+ { Name sc; NameOfCoef ic; Function BF_GroupOfEdges;
+ Support Dom_FW; Entity GroupsOfEdgesOf[Cut1H]; }
+
+ If(FE_Order == 2) // hierarchical basis functions for 2nd order interpolation
+
+ { Name se2; NameOfCoef he2; Function BF_Edge_2E;
+ Support Dom_FW; Entity EdgesOf[All, Not NonTerm]; }
+ // Reduced 2nd order
+ { Name se3a ; NameOfCoef he3a ; Function BF_Edge_3F_a ;
+ Support Dom_FW ; Entity FacetsOf[All, Not NonTerm] ; }
+ { Name se3b ; NameOfCoef he3b ; Function BF_Edge_3F_b ;
+ Support Dom_FW ; Entity FacetsOf[All, Not NonTerm] ; }
+ // 2nd order
+ { Name se4; NameOfCoef he4; Function BF_Edge_4E;
+ Support Dom_FW ; Entity EdgesOf[All, Not NonTerm]; }
+
+ { Name sn2; NameOfCoef phin2; Function BF_GradNode_2E;
+ Support Dom_FW ; Entity EdgesOf[BoundaryNotTerminal, Not SurfAirInf]; }
+
+ EndIf
+
+ Else
+ { Name se; NameOfCoef he; Function BF_Edge;
+ Support Dom_FW; Entity EdgesOf[All, Not BoundaryNotTerminal]; }
+
+ { Name sn; NameOfCoef phin; Function BF_GradNode;
+ Support Dom_FW;
+ //Entity NodesOf[BoundaryNotTerminal, Not Cut1H]; }
+ Entity NodesOf[BoundaryNotTerminal]; }
+
+ { Name sc; NameOfCoef ic; Function BF_GroupOfEdges;
+ Support Dom_FW; Entity GroupsOfEdgesOf[Cut1H]; }
+
+ If(FE_Order == 2) // hierarchical basis functions for 2nd order interpolation
+ // Edge BFs => FEorder = {0.5, Choices{0.5="Lowest",1="1st order",1.5="Reduced 2nd order",2="2nd order"}
+ // 1st order
+ { Name se2; NameOfCoef ee2; Function BF_Edge_2E;
+ Support Dom_FW; Entity EdgesOf[All, Not BoundaryNotTerminal]; }
+ // Reduced 2nd order
+ { Name se3a ; NameOfCoef ee3a ; Function BF_Edge_3F_a ;
+ Support Dom_FW ; Entity FacetsOf[All, Not BoundaryNotTerminal] ; }
+ { Name se3b ; NameOfCoef ee3b ; Function BF_Edge_3F_b ;
+ Support Dom_FW ; Entity FacetsOf[ All, Not BoundaryNotTerminal] ; }
+ // 2nd order
+ { Name se4; NameOfCoef ee4; Function BF_Edge_4E;
+ Support Dom_FW ; Entity EdgesOf[All, Not BoundaryNotTerminal]; }
+
+ // GradNode BFs
+ { Name sn2; NameOfCoef vn2; Function BF_GradNode_2E;
+ Support Dom_FW ; Entity EdgesOf[BoundaryNotTerminal]; }
+ EndIf
+ EndIf
+ }
+
+ GlobalQuantity {
+ { Name TerminalCurrent; Type AliasOf ; NameOfCoef ic; }
+ { Name TerminalPotential; Type AssociatedWith; NameOfCoef ic; }
+ }
+
+ Constraint {
+ { NameOfCoef TerminalPotential; EntityType GroupsOfEdgesOf;
+ NameOfConstraint SetTerminalPotentialH; }
+ { NameOfCoef TerminalCurrent; EntityType GroupsOfEdgesOf;
+ NameOfConstraint SetTerminalCurrentH; }
+ }
+ }
+}
+
+Function{
+ Omega = 2*Pi*Freq;
+ jOmega[] = Complex[0,1]*Omega;
+ jOmega2[] = -Omega^2;
+
+
+ tensSuma[] = sigma[]/jOmega[] + epsilon[];
+ alpha[] = 1/tensSuma[];
+
+ // alphaIm[] = jOmega[]*sigma[]/(sigma[]^2+Omega^2*epsilon[]^2);
+ // alphaRe[] = Omega^2*epsilon[]/(sigma[]^2+Omega^2*epsilon[]^2);
+ // alphaImSigmaZero[] = 0;
+ // alphaReSigmaZero[] = 1/epsilon[];
+
+}
+
+Formulation {
+
+ { Name FullWave_H_ece; Type FemEquation;
+ Quantity {
+ { Name h; Type Local; NameOfSpace Hcurl_H; }
+ { Name V; Type Global; NameOfSpace Hcurl_H[TerminalPotential]; }
+ { Name I; Type Global; NameOfSpace Hcurl_H[TerminalCurrent]; }
+ }
+ Equation {
+ // multiplying with j omega
+ // \int_D \alpha * curl(\vec{H}) \cdot curl(\vec{h}) dv
+
+ Galerkin { [ alpha[] * Dof{d h} , {d h} ]; In Vol_FW; Jacobian Vol; Integration Int; }
+
+ // \int_D j*\omega*(j*\omega*\mu) \vec{H} \cdot \vec{h} dv
+ Galerkin { [ jOmega2[] * mu[] * Dof{h} , {h} ]; In Vol_FW; Jacobian Vol; Integration Int; }
+
+ // sum (Vk ik); for k - all terminals so that the terminal voltages will be computed as well
+ GlobalTerm { [ jOmega[]*Dof{V}, {I} ]; In Cut1H; }
+
+ If (_use_radiantECE && Flag_SilverMuller)
+ Galerkin { DtDof [ Sqrt[1/epsilonNu[]] * ( Normal[] /\ Dof{h} ) /\ Normal[] , {h} ];
+ In SurfAirInf; Jacobian Sur; Integration Int; }
+ EndIf
+
+ }
+ }
+
+ { Name FullWave_H_ece_G; Type FemEquation;
+ Quantity {
+ { Name h; Type Local; NameOfSpace Hcurl_H; }
+ { Name V; Type Global; NameOfSpace Hcurl_H[TerminalPotential]; }
+ { Name I; Type Global; NameOfSpace Hcurl_H[TerminalCurrent]; }
+ }
+ Equation {
+ // without multiplying with j omega
+ Galerkin { [ (1/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[])) * Dof{d h} , {d h} ]; In Vol_FW; Jacobian Vol; Integration Int; }
+ Galerkin { [ Complex[0,1]*2*Pi*Freq* mu[] * Dof{h} , {h} ]; In Vol_FW; Jacobian Vol; Integration Int; }
+
+ // sum (Vk ik); for k - all terminals so that the terminal voltages will be computed as well
+ GlobalTerm { [ Dof{V}, {I} ]; In Cut1H; }
+ }
+ }
+}