From 08403f8d9039eb05fc5e4bbe721eae1ed8d8807d Mon Sep 17 00:00:00 2001
From: "Ruth Sabariego (U0089683)" <Ruth.Sabariego@esat.kuleuven.be>
Date: Mon, 13 Mar 2023 15:43:19 +0100
Subject: [PATCH] Fullwave formulations with ECE constraints. Two academic test
 cases: a cylinder and a coaxial cable

---
 cylinderCoax/cylinderCoax.geo                 |  17 ++
 cylinderCoax/cylinderCoax.pro                 |  17 ++
 cylinderCoax/cylinderCoax_data.pro            | 120 +++++++++
 ...lWave_E_ece_MIMO2terminals_secondOrder.pro |  92 +++++++
 .../FullWave_E_ece_SISO_secondOrder.pro       |  89 ++++++
 ...2terminals_cplx_2D_and_AXI_secondOrder.pro |  87 ++++++
 ...ece_MIMO2terminals_cplx_3D_secondOrder.pro |  87 ++++++
 ...H_ece_SISO_cplx_2D_and_AXI_secondOrder.pro | 104 +++++++
 ...ullWave_H_ece_SISO_cplx_3D_secondOrder.pro |  61 +++++
 .../formulations/Jacobian_and_Integration.pro |  68 +++++
 ...Formulation_FullWave_E_ece_secondOrder.pro |  80 ++++++
 ...Wave_H_ece_cplx_2D_and_AXI_secondOrder.pro |  72 +++++
 ...ion_FullWave_H_ece_cplx_3D_secondOrder.pro | 116 ++++++++
 cylinderCoax/geo_pro/coax2Daxi.geo            |  55 ++++
 cylinderCoax/geo_pro/coax2Daxi.pro            | 253 ++++++++++++++++++
 cylinderCoax/geo_pro/coax2Daxi_data.pro       | 178 ++++++++++++
 cylinderCoax/geo_pro/coax3D.geo               |  99 +++++++
 cylinderCoax/geo_pro/coax3D.pro               | 217 +++++++++++++++
 cylinderCoax/geo_pro/coax3D_data.pro          | 174 ++++++++++++
 cylinderCoax/geo_pro/cylinder2Daxi.geo        | 145 ++++++++++
 cylinderCoax/geo_pro/cylinder2Daxi.pro        | 175 ++++++++++++
 cylinderCoax/geo_pro/cylinder2Daxi_data.pro   | 202 ++++++++++++++
 cylinderCoax/geo_pro/cylinder3D.geo           |  89 ++++++
 cylinderCoax/geo_pro/cylinder3D.pro           | 160 +++++++++++
 cylinderCoax/geo_pro/cylinder3D_data.pro      | 167 ++++++++++++
 25 files changed, 2924 insertions(+)
 create mode 100644 cylinderCoax/cylinderCoax.geo
 create mode 100644 cylinderCoax/cylinderCoax.pro
 create mode 100644 cylinderCoax/cylinderCoax_data.pro
 create mode 100644 cylinderCoax/formulations/FullWave_E_ece_MIMO2terminals_secondOrder.pro
 create mode 100644 cylinderCoax/formulations/FullWave_E_ece_SISO_secondOrder.pro
 create mode 100644 cylinderCoax/formulations/FullWave_H_ece_MIMO2terminals_cplx_2D_and_AXI_secondOrder.pro
 create mode 100644 cylinderCoax/formulations/FullWave_H_ece_MIMO2terminals_cplx_3D_secondOrder.pro
 create mode 100644 cylinderCoax/formulations/FullWave_H_ece_SISO_cplx_2D_and_AXI_secondOrder.pro
 create mode 100644 cylinderCoax/formulations/FullWave_H_ece_SISO_cplx_3D_secondOrder.pro
 create mode 100644 cylinderCoax/formulations/Jacobian_and_Integration.pro
 create mode 100644 cylinderCoax/formulations/only_FunctionSpace_and_Formulation_FullWave_E_ece_secondOrder.pro
 create mode 100644 cylinderCoax/formulations/only_FunctionSpace_and_Formulation_FullWave_H_ece_cplx_2D_and_AXI_secondOrder.pro
 create mode 100644 cylinderCoax/formulations/only_FunctionSpace_and_Formulation_FullWave_H_ece_cplx_3D_secondOrder.pro
 create mode 100644 cylinderCoax/geo_pro/coax2Daxi.geo
 create mode 100644 cylinderCoax/geo_pro/coax2Daxi.pro
 create mode 100644 cylinderCoax/geo_pro/coax2Daxi_data.pro
 create mode 100644 cylinderCoax/geo_pro/coax3D.geo
 create mode 100644 cylinderCoax/geo_pro/coax3D.pro
 create mode 100644 cylinderCoax/geo_pro/coax3D_data.pro
 create mode 100644 cylinderCoax/geo_pro/cylinder2Daxi.geo
 create mode 100644 cylinderCoax/geo_pro/cylinder2Daxi.pro
 create mode 100644 cylinderCoax/geo_pro/cylinder2Daxi_data.pro
 create mode 100644 cylinderCoax/geo_pro/cylinder3D.geo
 create mode 100644 cylinderCoax/geo_pro/cylinder3D.pro
 create mode 100644 cylinderCoax/geo_pro/cylinder3D_data.pro

diff --git a/cylinderCoax/cylinderCoax.geo b/cylinderCoax/cylinderCoax.geo
new file mode 100644
index 0000000..46b494a
--- /dev/null
+++ b/cylinderCoax/cylinderCoax.geo
@@ -0,0 +1,17 @@
+Include "cylinderCoax_data.pro";
+
+If (Flag_model == 1)     // cylinder
+  If (Flag_3Dmodel == 0) // 2D axi
+    Include "./geo_pro/cylinder2Daxi.geo";
+  Else // 3D
+    Include "./geo_pro/cylinder3D.geo";
+  EndIf
+EndIf
+
+If (Flag_model == 2)  // coaxial cable
+  If (Flag_3Dmodel == 0) // 2D axi
+    Include "./geo_pro/coax2Daxi.geo";
+  Else // 3D
+    Include "./geo_pro/coax3D.geo";
+  EndIf
+EndIf
diff --git a/cylinderCoax/cylinderCoax.pro b/cylinderCoax/cylinderCoax.pro
new file mode 100644
index 0000000..7bacc60
--- /dev/null
+++ b/cylinderCoax/cylinderCoax.pro
@@ -0,0 +1,17 @@
+Include "cylinderCoax_data.pro";
+
+If (Flag_model == 1)     // cylinder
+  If (Flag_3Dmodel == 0) // 2D axi
+    Include "./geo_pro/cylinder2Daxi.pro";
+  Else // 3D
+    Include "./geo_pro/cylinder3D.pro";
+  EndIf
+EndIf
+
+If (Flag_model == 2)  // coaxial cable
+  If (Flag_3Dmodel == 0) // 2D axi
+    Include "./geo_pro/coax2Daxi.pro";
+  Else // 3D
+    Include "./geo_pro/coax3D.pro";
+  EndIf
+EndIf
diff --git a/cylinderCoax/cylinderCoax_data.pro b/cylinderCoax/cylinderCoax_data.pro
new file mode 100644
index 0000000..28276da
--- /dev/null
+++ b/cylinderCoax/cylinderCoax_data.pro
@@ -0,0 +1,120 @@
+// Cylinder and Coaxial Cable with ECE
+// uncomment only one of the lines TEST = ..... below
+
+// The geometries are in separate .geo files 2.5D / 3D   cylinder/coax => 4 geo files in the geo folder
+//                                           
+// They use separate pro files => see pro folder
+
+// The following tests should work without errors
+// TEST code of the test 7 digits - their meaning is given in the following order:
+// 1234567
+// 1 model ID: 1 for cylinder, 2 for coaxial cable
+// 2 is the submodel number 
+// 3 is 0 for voltage excitation and 1 for current excitation of terminal no 1 (cylinder - there is only one excited terminal)
+// 4 is 0 for voltage excitation and 1 for current excitation of terminal no 1 (it has a meaning only for the coaxial cable example)
+// 5 is for the model flag 0 for 2.5D and 1 for 3D
+// 6 is for the formulation 0 for E and 1 for H
+// 7 id for FE order 1 or  2
+
+DefineConstant [ // allowing an external call (-set_number), e.g. from Matlab, python or a shell	
+    TEST = 0  // for no predefined tests
+	
+	// TEST code of the test 7 digits - their meaning is given in the following order:
+	// 1234567
+	// 1 - model ID: 1 for cylindrical, 2 for coax
+	// 2 - is the submodel ID
+	// 3 is 0 for voltage excitation and 1 for current excitation - for the first excited terminal (the only one for the cylinder test)
+	// 4 is 0 for voltage excitation and 1 for current excitation - for the second excited terminal (relevant only for the coax test)
+	// 5 is for the model flag 0 for 2.5D and 1 for 3D
+	// 6 is for the formulation 0 for E and 1 for H
+	// 7 id for FE order 1 or  2
+
+	// TESTS - you can look in the s1p files in the res folder
+	//------ cylinder HF, v1, 2.5D ------
+	//TEST = 1100001    // 1) cylinder, v1, voltage excitation, - , 2.5D, formulation in E, 1st order elements  
+	//TEST = 1100002    // 2) cylinder, v1, voltage excitation, - , 2.5D, formulation in E, 2nd order elements 
+	//TEST = 1100011    // 3) cylinder, v1, voltage excitation, - , 2.5D, formulation in H, 1st order elements 
+	//TEST = 1100012    // 4) cylinder, v1, voltage excitation, - , 2.5D, formulation in H, 2nd order elements 
+    //TEST = 1110001    // 5) cylinder, v1, current excitation, - , 2.5D, formulation in E, 1st order elements  
+    //TEST = 1110002    // 6) cylinder, v1, current excitation, - , 2.5D, formulation in E, 2nd order elements
+	//TEST = 1110011    // 7) cylinder, v1, current excitation, - , 2.5D, formulation in H, 1st order elements   
+	//TEST = 1110012    // 8) cylinder, v1, current excitation, - , 2.5D, formulation in H, 2nd order elements 
+    //------ coax HF, 2.5D ------
+	//TEST = 2100001    // 9) coax HF, voltage excitation, - , 2.5D, formulation in E, 1st order elements       // needed for fig 5 in the SCEE 2022 paper
+	//TEST = 2100002    // 10) coax HF, voltage excitation, - , 2.5D, formulation in E, 2nd order elements      // needed for fig 5 in the SCEE 2022 paper
+	//TEST = 2100011    // 11) coax HF, voltage excitation, - , 2.5D, formulation in H, 1st order elements      // needed for fig 5 in the SCEE 2022 paper
+	//TEST = 2100012    // 12) coax HF, voltage excitation, - , 2.5D, formulation in H, 2nd order elements      // needed for fig 5 in the SCEE 2022 paper
+    //TEST = 2111001    // 13) coax HF, current excitation, - , 2.5D, formulation in E, 1st order elements  
+    //TEST = 2111002    // 14) coax HF, current excitation, - , 2.5D, formulation in E, 2nd order elements
+	//TEST = 2111011    // 15) coax HF, current excitation, - , 2.5D, formulation in H, 1st order elements   
+	//TEST = 2111012    // 16) coax HF, current excitation, - , 2.5D, formulation in H, 2nd order elements 
+    //------ coax LF, 2.5D ------
+	//TEST = 2200001    // 17) coax LF, voltage excitation, - , 2.5D, formulation in E, 1st order elements      // needed for fig 6 in the SCEE 2022 paper
+	//TEST = 2200002    // 18) coax LF, voltage excitation, - , 2.5D, formulation in E, 2nd order elements      // needed for fig 6 in the SCEE 2022 paper
+	//TEST = 2200011    // 19) coax LF, voltage excitation, - , 2.5D, formulation in H, 1st order elements      // needed for fig 6 in the SCEE 2022 paper
+	//TEST = 2200012    // 20) coax LF, voltage excitation, - , 2.5D, formulation in H, 2nd order elements      // needed for fig 6 in the SCEE 2022 paper
+    //TEST = 2211001    // 21) coax LF, current excitation, - , 2.5D, formulation in E, 1st order elements  
+    //TEST = 2211002    // 22) coax LF, current excitation, - , 2.5D, formulation in E, 2nd order elements
+	//TEST = 2211011    // 23) coax LF, current excitation, - , 2.5D, formulation in H, 1st order elements   
+	//TEST = 2211012    // 24) coax LF, current excitation, - , 2.5D, formulation in H, 2nd order elements 
+
+	//------ 3D
+	//------ cylinder, v1, 3D ------
+	//TEST = 1100101    // 101) cylinder, v1, voltage excitation, - , 3D, formulation in E, 1st order elements  
+	//TEST = 1100102    // 102) cylinder, v1, voltage excitation, - , 3D, formulation in E, 2nd order elements 
+	//TEST = 1100111    // 103) cylinder, v1, voltage excitation, - , 3D, formulation in H, 1st order elements 
+	//TEST = 1100112    // 104) cylinder, v1, voltage excitation, - , 3D, formulation in H, 2nd order elements 
+    //TEST = 1110101    // 105) cylinder, v1, current excitation, - , 3D, formulation in E, 1st order elements   // needed for fig 1 in the SCEE 2022 paper
+    //TEST = 1110102    // 106) cylinder, v1, current excitation, - , 3D, formulation in E, 2nd order elements
+	//TEST = 1110111    // 107) cylinder, v1, current excitation, - , 3D, formulation in H, 1st order elements   // needed for fig 1 in the SCEE 2022 paper
+	//TEST = 1110112    // 108) cylinder, v1, current excitation, - , 3D, formulation in H, 2nd order elements 
+];
+
+
+
+If (TEST == 0)  // some default setting if no predefined test is used
+	MODEL_ID = 0;   // no model chosen
+	FLAG_MODEL = 0; // 2D axi
+	
+	TERMINAL_EXC = 1;
+	TERMINAL_EXC2 = 1; 
+	FORMULATION = 0;
+	FE_ORDER = 1;
+	
+	MODEL_ID = 1 ;               // 1 for cylindrical, 2 for coax
+	MODEL_ID2 =1 ;              // 1 vor version v1 - not really used for the moment
+	TERMINAL_EXC = 0;        // 0 for voltage excitation, 1 for current excitation
+	TERMINAL_EXC2 = 0;       // second terminal - relevant only for the coax test
+	FLAG_MODEL = 1;          // 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
+EndIf
+
+/* 2.5D */
+// Include "./tests/tests_1to8_data.pro";         /* cylinder - version v1 (SCEE 2020, JMI) - 2.5D */
+// Include "./tests/tests_9to16_data.pro";        /* coaxial cable [Wu04] with conductor inside, HF (SCEE 2022)  - 2.5D */
+// Include "./tests/tests_17to24_data.pro";       /* coaxial cable [Wu04] with conductor inside, LF (SCEE 2022)  - 2.5D */
+
+/* 3D */
+// Include "./tests/tests_101to108_data.pro";    /* cylinder - version v1 (SCEE 2020, 2022, JMI) - 3D */
+//Include "./tests/tests_109to116_data.pro";       /* coaxial cable [Wu04] with conductor inside - 3D */
+
+// If TEST == 0 , all the tests above were commented 
+/* _______Choice 1_______  chose the model 2.5D (2D axi) or 3D  */
+DefineConstant[
+	Flag_3Dmodel = {0, Choices {0,1}, 
+		Name "{Model/1Use 3D model?", Highlight "Pink"}
+	ang_section_degrees = 360
+	Flag_model = {1, Choices{1="Cylinder (SISO)",2="Coaxial cable (MIMO)"}, 
+		Name "{Model/2Choose model", Highlight "Tomato"}
+];	
+
+modelDim = (!Flag_3Dmodel)? 2: 3;
+Flag_Axi = (!Flag_3Dmodel)? 1:-1; // 2.5 D; if this flag is zero, this would mean a 2D problem, not an axisymmetric one
+h2Ddepth = (!Flag_3Dmodel)? 2*Pi:360/ang_section_degrees;
+
+
+// just to be sure that the parameters are correct
+Printf("========================= MAIN PARAMETERS =========================");
+Printf("*** *** *** ==> Flag_model   (1-Cylinder, 2-Coaxial cable, 0-Not selected) = %e",Flag_model);
+Printf("*** *** *** ==> Flag_3Dmodel (0-for 2.5D, 1 for 3D) = %e",Flag_3Dmodel);
diff --git a/cylinderCoax/formulations/FullWave_E_ece_MIMO2terminals_secondOrder.pro b/cylinderCoax/formulations/FullWave_E_ece_MIMO2terminals_secondOrder.pro
new file mode 100644
index 0000000..bdb3434
--- /dev/null
+++ b/cylinderCoax/formulations/FullWave_E_ece_MIMO2terminals_secondOrder.pro
@@ -0,0 +1,92 @@
+/* FullWave_E_ece_MIMO2terminals.pro  (old name "formulationFW_E_eceBC_v6_MIMO2terminals")  
+
+ Problem independent pro file
+ for passive, linear domains in Full Wave with 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) and the other one is set to zero.
+ Frequency domain simulations.
+*/
+
+Include "Jacobian_and_Integration.pro"
+Include "only_FunctionSpace_and_Formulation_FullWave_E_ece_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
+      Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
+    }
+  }
+
+}
+
+PostProcessing {
+
+  { Name FW_E_ece; NameOfFormulation FullWave_E_ece;
+    Quantity {
+      { Name E; Value {
+            Local { [ {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 Vterminal; Value {
+        Term { [ -{V} ]; In Terminal1; } 
+        Term { [ -{V} ]; In Terminal2; } 
+        }
+      }
+
+      { Name Iterminal; Value { 
+        Term { [ -1*{I}*h2Ddepth ]; In Terminal1; } 
+        Term { [ -1*{I}*h2Ddepth ]; In Terminal2; } 
+        }  
+      } 
+      
+      { Name VnotTerminals; Value {
+          Local { [ -{dInv dv} ]; In Vol_FW; Jacobian Vol; }
+        }
+      }
+	  
+      { Name B; Value {
+          Local { [ CompZ[{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 Hz; Value {
+        Local { [ nu[]*CompZ[{d e}]/(2*Pi*Freq*Complex[0,1]) ]; In Vol_FW; Jacobian Vol; }
+        }
+      }
+    }
+  }
+}
diff --git a/cylinderCoax/formulations/FullWave_E_ece_SISO_secondOrder.pro b/cylinderCoax/formulations/FullWave_E_ece_SISO_secondOrder.pro
new file mode 100644
index 0000000..d615124
--- /dev/null
+++ b/cylinderCoax/formulations/FullWave_E_ece_SISO_secondOrder.pro
@@ -0,0 +1,89 @@
+/* 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_ece_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
+			Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
+		}
+	}
+}
+
+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 Vterminal; Value {
+				  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 {
+        Term { [ -{I}*h2Ddepth ]; In Sur_Terminals_FWece; }
+				}
+			}  // h2Ddepth is the depth for 2D problems,  and it has to be 1 for 3D problems
+
+			// 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; }
+				}
+			}
+			
+		}
+	}
+ }
diff --git a/cylinderCoax/formulations/FullWave_H_ece_MIMO2terminals_cplx_2D_and_AXI_secondOrder.pro b/cylinderCoax/formulations/FullWave_H_ece_MIMO2terminals_cplx_2D_and_AXI_secondOrder.pro
new file mode 100644
index 0000000..16badf7
--- /dev/null
+++ b/cylinderCoax/formulations/FullWave_H_ece_MIMO2terminals_cplx_2D_and_AXI_secondOrder.pro
@@ -0,0 +1,87 @@
+/* FullWave_H_ece_SISO_cplx_3D.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_ece_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
+    Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
+    }
+  }
+
+}
+
+PostProcessing {
+ { Name FW_H_ece; NameOfFormulation FullWave_H_ece;
+    Quantity {
+      { Name H; Value {
+        Local { [ {h} ]; In Dom_FW; Jacobian Vol; }
+        }
+      }
+
+      { Name Hz ; Value { 
+        Local{ [ CompZ[{h}] ] ; In Dom_FW ; Jacobian Vol ; }
+        }
+      }
+
+      { Name E; Value {
+        Local { [ {d h}/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[]) ]; In Dom_FW; Jacobian Vol; }
+        }
+      }
+      
+      { Name Ey; Value {
+          Local { [ CompY[{d h}/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[])] ]; In Dom_FW; Jacobian Vol; }
+        }
+      }
+
+      { Name J; Value {
+        Local { [ {d h}*sigma[]/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[]) ]; In Dom_FW; Jacobian Vol; }
+        }
+      }
+    
+      { Name rmsJ; Value {
+        Local { [ Norm[{d h}*sigma[]/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[])] ]; In Dom_FW; Jacobian Vol; }
+        }
+      }
+      
+      { Name rmsH; Value {
+        Local { [ Norm[{h}] ]; In Dom_FW; Jacobian Vol; }
+        }
+      }      
+
+      { Name Vterminal; Value {
+        Term { [  {V} ]; In Cut1Hterminal1; }  
+        Term { [ -{V} ]; In Cut1Hterminal2; }
+        }
+      }
+      
+      { Name Iterminal; Value { 
+        Term { [ -{I}*h2Ddepth ]; In Cut1Hterminal1; } 
+        Term { [  {I}*h2Ddepth ]; In Cut1Hterminal2; } 
+        }     
+      }
+    } 
+  }
+}
diff --git a/cylinderCoax/formulations/FullWave_H_ece_MIMO2terminals_cplx_3D_secondOrder.pro b/cylinderCoax/formulations/FullWave_H_ece_MIMO2terminals_cplx_3D_secondOrder.pro
new file mode 100644
index 0000000..b3ee31c
--- /dev/null
+++ b/cylinderCoax/formulations/FullWave_H_ece_MIMO2terminals_cplx_3D_secondOrder.pro
@@ -0,0 +1,87 @@
+/* FullWave_H_ece_SISO_cplx_3D.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_ece_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];              
+      Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
+    }
+  }
+
+}
+
+PostProcessing {
+  { Name FW_H_ece; NameOfFormulation FullWave_H_ece;
+    Quantity {
+      { Name H; Value {
+          Local { [ {h} ]; In Dom_FW; Jacobian Vol; } 
+        }
+      }
+
+      { Name Hz ; Value { 
+		      Local{ [ CompZ[{h}] ] ; In Dom_FW ; Jacobian Vol ; }
+        }
+      }
+
+	    { Name E; Value {
+          Local { [ {d h}/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[]) ]; In Dom_FW; Jacobian Vol; }
+        }
+      }
+
+	    { Name J; Value {
+          Local { [ {d h}*sigma[]/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[]) ]; In Dom_FW; Jacobian Vol; }
+        }
+      }
+
+	    { Name rmsJ; Value {
+          Local { [ Norm[{d h}*sigma[]/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[])] ]; In Dom_FW; Jacobian Vol; }
+        }
+      }
+	  
+	    { Name rmsH; Value {
+          Local { [ Norm[{h}] ]; In Dom_FW; Jacobian Vol; }
+        }
+      }
+
+	    { Name Ey; Value {
+          Local { [ CompY[{d h}/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[])] ]; In Dom_FW; Jacobian Vol; }
+        }
+      }
+
+	    { Name Vterminal; Value {
+		    Term { [  {V} ]; In Cut1Hterminal1; }
+		    Term { [ -{V} ]; In Cut1Hterminal2; }
+        }
+      }
+	  
+      { Name Iterminal; Value { 
+          Term { [ -{I}*h2Ddepth ]; In Cut1Hterminal1; }   
+          Term { [ -{I}*h2Ddepth ]; In Cut1Hterminal2; }   
+          // h is the depth for 2D problems,  and it has to be 1 for 3D problems
+  	    }
+      }    
+    }
+  } 
+}
diff --git a/cylinderCoax/formulations/FullWave_H_ece_SISO_cplx_2D_and_AXI_secondOrder.pro b/cylinderCoax/formulations/FullWave_H_ece_SISO_cplx_2D_and_AXI_secondOrder.pro
new file mode 100644
index 0000000..4a776b5
--- /dev/null
+++ b/cylinderCoax/formulations/FullWave_H_ece_SISO_cplx_2D_and_AXI_secondOrder.pro
@@ -0,0 +1,104 @@
+/* FullWave_H_ece_SISO_cplx_2D.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_ece_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];               
+			Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
+		}
+	}
+}
+
+PostProcessing {
+
+	{ Name FW_H_ece; NameOfFormulation FullWave_H_ece;
+		Quantity {
+			{ Name H;    Value { Local { [ {h} ];       In Vol_FW; Jacobian Vol; }}} 
+			{ Name rmsH; Value { Local { [ Norm[{h}] ]; In Vol_FW; Jacobian Vol; }}}
+			{ Name Hz ;  Value { Local{ [ CompZ[{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 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 Ey;	Value {	Local { [ CompY[{d h}/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[])] ]; In Vol_FW; Jacobian Vol; }}}
+
+			{ Name Vterminal;	Value {	Term { [ -1*{V} ]; In Cut1H; }}}
+
+			{ Name Rin;	Value {Term { [ Re[-1*{V}] ]; In Cut1H; }}}
+			{ Name Xin;	Value {Term { [ Im[-1*{V}] ]; In Cut1H; }}}
+			{ Name AbsZin; Value { 
+          Term { [ Norm[-1*{V}] ]; In Cut1H; } 
+				}
+			}
+			{ Name ArgZin; Value { 
+          Term { [Atan2[Im[-1*{V} ],Re[ -1*{V} ]]]; In Cut1H; } 
+				}
+			}
+	  
+			{ Name I;         	     
+				Value { 
+					If(Flag_Axi == 0) // 2D 
+						Term { [ -2*{I}*h2Ddepth ]; In Cut1H; } 
+					Else  // 2D axisymmetric
+						Term { [ -{I}*h2Ddepth ]; In Cut1H; } 
+					EndIf
+				}   // h is the depth for 2D problems,  and it has to be 1 for 3D problems
+			} 
+
+			{ 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 { 
+          Term { [Atan2[Im[ -{I}*h2Ddepth ],Re[ -{I}*h2Ddepth ]]]; In Cut1H; } 
+				}
+			} 	    
+		}
+	}
+}  
+ 
+ 
+	  
+	
diff --git a/cylinderCoax/formulations/FullWave_H_ece_SISO_cplx_3D_secondOrder.pro b/cylinderCoax/formulations/FullWave_H_ece_SISO_cplx_3D_secondOrder.pro
new file mode 100644
index 0000000..e53aafd
--- /dev/null
+++ b/cylinderCoax/formulations/FullWave_H_ece_SISO_cplx_3D_secondOrder.pro
@@ -0,0 +1,61 @@
+/* FullWave_H_ece_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_ece_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]; 
+      Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
+    }
+  }
+}
+
+
+PostProcessing {
+
+  { Name FW_H_ece; NameOfFormulation FullWave_H_ece;
+    Quantity {
+      { Name H;    Value { Local { [ {h} ];       In Vol_FW; Jacobian Vol; }}} 
+      { Name rmsH; Value { Local { [ Norm[{h}] ]; In Vol_FW; Jacobian Vol; }}}
+      { Name Hz ;  Value { Local{ [ CompZ[{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 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 Ey;	Value {	Local { [ CompY[{d h}/(sigma[]+Complex[0,1]*2*Pi*Freq*epsilon[])] ]; In Vol_FW; Jacobian Vol; }}}
+
+      { Name Vterminals;	Value {	Term { [ -1*{V} ]; In Cut1H; }}}
+
+      { Name Rin;	Value {Term { [ Re[-1*{V}] ]; In Cut1H; }}}
+      { Name Xin;	Value {Term { [ Im[-1*{V}] ]; In Cut1H; }}}
+      { Name AbsZin;   	Value { Term { [Sqrt[Im[ -1*{V} ]^2 + Re[ -1*{V} ]^2]]; In Cut1H; }}}
+      { Name ArgZin; 		Value { Term { [Atan2[Im[ -1*{V} ],Re[ -1*{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/cylinderCoax/formulations/Jacobian_and_Integration.pro b/cylinderCoax/formulations/Jacobian_and_Integration.pro
new file mode 100644
index 0000000..42eb12f
--- /dev/null
+++ b/cylinderCoax/formulations/Jacobian_and_Integration.pro
@@ -0,0 +1,68 @@
+Jacobian {
+	{ Name Vol;
+		Case {
+			If(Flag_Axi) 
+				{ 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/cylinderCoax/formulations/only_FunctionSpace_and_Formulation_FullWave_E_ece_secondOrder.pro b/cylinderCoax/formulations/only_FunctionSpace_and_Formulation_FullWave_E_ece_secondOrder.pro
new file mode 100644
index 0000000..88d5adb
--- /dev/null
+++ b/cylinderCoax/formulations/only_FunctionSpace_and_Formulation_FullWave_E_ece_secondOrder.pro
@@ -0,0 +1,80 @@
+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
+
+				{ 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		 
+		 
+	  
+			{ Name sf; NameOfCoef vf; Function BF_GradGroupOfNodes; 
+				Support Dom_FW ; Entity GroupsOfNodesOf[Sur_Terminals_FWece]; }
+			 
+		}
+		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; }
+	  
+		}
+	}
+}
+ 
+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; }
+
+			//    alternative // j*\omega*sum (vk Ik);  for k - current excited terminals
+			//    GlobalTerm {DtDof [ -Dof{I} , {V} ]; In SurBCec; }
+
+			// 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; }  
+		}
+	}
+}
diff --git a/cylinderCoax/formulations/only_FunctionSpace_and_Formulation_FullWave_H_ece_cplx_2D_and_AXI_secondOrder.pro b/cylinderCoax/formulations/only_FunctionSpace_and_Formulation_FullWave_H_ece_cplx_2D_and_AXI_secondOrder.pro
new file mode 100644
index 0000000..a6606a1
--- /dev/null
+++ b/cylinderCoax/formulations/only_FunctionSpace_and_Formulation_FullWave_H_ece_cplx_2D_and_AXI_secondOrder.pro
@@ -0,0 +1,72 @@
+//=========================================================
+// Function spaces
+//=========================================================
+
+FunctionSpace {
+  
+	{ Name Hcurl_H; Type Form1P;  // Hspace for 2D problems with H perpendicular to the plane
+		BasisFunction {
+      { 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		  
+		}
+
+		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
+		}
+	}
+}
+
+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 { [ alphaRe[] * Dof{d h} , {d h} ];      In Vol_FW; Jacobian Vol; Integration Int; }
+			Galerkin { [ alphaIm[] * 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; }
+		}
+	}
+}
diff --git a/cylinderCoax/formulations/only_FunctionSpace_and_Formulation_FullWave_H_ece_cplx_3D_secondOrder.pro b/cylinderCoax/formulations/only_FunctionSpace_and_Formulation_FullWave_H_ece_cplx_3D_secondOrder.pro
new file mode 100644
index 0000000..4478882
--- /dev/null
+++ b/cylinderCoax/formulations/only_FunctionSpace_and_Formulation_FullWave_H_ece_cplx_3D_secondOrder.pro
@@ -0,0 +1,116 @@
+//=========================================================
+// Function spaces
+//=========================================================
+
+FunctionSpace {
+
+	{ Name Hcurl_H; Type Form1;   // Hspace for 3D
+		BasisFunction {
+	
+				{ 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
+				
+    }
+	
+	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;
+
+  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 { [ alphaRe[] * Dof{d h} , {d h} ];      In Vol_FW; Jacobian Vol; Integration Int; }
+      Galerkin { [ alphaIm[] * Dof{d h} , {d h} ];      In Vol_FW; Jacobian Vol; Integration Int; }
+	      	  
+		/* 
+	  Galerkin { [ alphaRe[] * Dof{d h} , {d h} ];      In Vol_FW_SigmaNonZero; Jacobian Vol; Integration Int; }
+      Galerkin { [ alphaIm[] * Dof{d h} , {d h} ];      In Vol_FW_SigmaNonZero; Jacobian Vol; Integration Int; }
+      Galerkin { [ alphaReSigmaZero[] * Dof{d h} , {d h} ];  In Vol_FW_SigmaZero; Jacobian Vol; Integration Int; }
+      //Galerkin { [ alphaImSigmaZero[] * Dof{d h} , {d h} ];      In Vol_FW_SigmaZero; 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; }
+	  
+    }
+  }
+  
+  { 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; }
+    }
+  }
+}
diff --git a/cylinderCoax/geo_pro/coax2Daxi.geo b/cylinderCoax/geo_pro/coax2Daxi.geo
new file mode 100644
index 0000000..c829d7d
--- /dev/null
+++ b/cylinderCoax/geo_pro/coax2Daxi.geo
@@ -0,0 +1,55 @@
+Include "coax2Daxi_data.pro";
+
+// Definition of parameters for local mesh dimensions
+p0 = s*l/10; // characteristic length of mesh element
+
+// Definition of gemetrical points
+Point(1) = { 0, 0, 0, p0} ;
+Point(2) = { a, 0, 0, p0} ;
+Point(3) = { a, l, 0, p0} ;
+Point(4) = { 0, l, 0, p0} ;
+Point(5) = { b, 0, 0, p0} ;
+Point(6) = { b, l, 0, p0} ;
+
+// Definition of gemetrical lines
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {4,1};
+Line(5) = {2,5};
+Line(6) = {5,6};
+Line(7) = {3,6};
+
+// Definition of geometrical surfaces
+Line Loop(5) = {1, 2, 3, 4};
+Plane Surface(6) = {5};
+Line Loop(7) = {5, 6, -7, -2};
+Plane Surface(8) = {7};
+
+If(_use_transfinite)
+  Transfinite Line {2,4,6} = 9*s; // 8  discretizations
+  Transfinite Line {1,3,5,7} = 4*s; 
+  Transfinite Surface {6,8};
+  If (_use_recombine)
+    Recombine Surface{6,8};
+  EndIf
+EndIf
+
+/* Definition of Physical entities (surfaces, lines). The Physical
+   entities tell GMSH the elements and their associated region numbers
+   to save in the file 'Ishape2d.msh'. */
+
+Physical Surface ("Conductor", 100) = {6} ;   
+Physical Surface ("Air", 101) = {8} ;    
+
+Physical Line ("Terminal1", 121) = {1} ;          
+Physical Line ("Terminal2", 122) = {3} ;          
+Physical Line ("Ground",   120) = {6} ;          
+
+Physical Line ("LeftBoundary", 132) = {4} ;      
+Physical Line ("BottomBoundary", 133) = {5} ; 
+Physical Line ("TopBoundary", 134) = {7} ; 
+
+// For aestetics
+Recursive Color SkyBlue { Surface{8};}
+Recursive Color Red  { Surface{6};}
\ No newline at end of file
diff --git a/cylinderCoax/geo_pro/coax2Daxi.pro b/cylinderCoax/geo_pro/coax2Daxi.pro
new file mode 100644
index 0000000..6ac9de3
--- /dev/null
+++ b/cylinderCoax/geo_pro/coax2Daxi.pro
@@ -0,0 +1,253 @@
+Include "coax2Daxi_data.pro";
+
+Group{
+  Conductor = Region[100]; 
+  Air = Region[101];
+  
+  Ground = Region[120];    /* Boundary - various parts */
+  Terminal1 = Region[121];
+  Terminal2 = Region[122];
+  
+
+  BottomBoundary = Region[133];
+  TopBoundary = Region[134];
+  BoundaryNotTerminal = Region[{BottomBoundary,TopBoundary}];
+
+  If(Flag_Axi == 0) // 2D
+    LeftBoundary = Region[132];
+    Axis = Region[{}];
+    BoundaryNotTerminal += Region[{LeftBoundary}];
+  Else // axi
+    LeftBoundary = Region[{}];
+    Axis = Region[132];
+  EndIf
+  
+  Cut1Hterminal1 = Region[{BottomBoundary}];
+  Cut1Hterminal2 = Region[{TopBoundary}];
+  
+  Cut1H = Region[{Cut1Hterminal1,Cut1Hterminal2}];
+  
+  Sur_Terminals_FWece = Region[{Ground, Terminal1, Terminal2}]; // all terminals
+  
+  Vol_FW = Region[{Conductor,Air}]; // domain without the boundary
+  Sur_FW = Region[{Sur_Terminals_FWece, TopBoundary, BottomBoundary, LeftBoundary, Axis}];
+  Dom_FW = Region[ {Vol_FW, Sur_FW} ];
+
+  If (sigmaExt == 0) // used only with H formulation
+    Vol_FW_SigmaZero = Region[{Air}];
+    Vol_FW_SigmaNonZero = Region[{Conductor}];
+  Else
+    Vol_FW_SigmaZero = Region[{}];
+    Vol_FW_SigmaNonZero = Region[{Air,Conductor}];
+  EndIf
+  
+}
+
+Function{
+   /* Material properties are defined piecewise, for elementary groups */
+   /* the values in the right hand side are defined in Ishape2d_data.pro */
+   
+   nu[#{Conductor,Terminal1,Terminal2,Axis}] = nuInt;  // magnetic reluctivity [m/H]
+   mu[#{Conductor,Terminal1,Terminal2,Axis}] = muInt;  // permeability [H/m]
+   epsilon[#{Conductor,Terminal1,Terminal2,Axis}] = epsInt;       // permittivity [F/m]
+   sigma[#{Conductor,Terminal1,Terminal2,Axis}]  = sigmaInt;      // conductivity [S/m]
+   
+   nu[#{Air,BoundaryNotTerminal}] = nuExt;  // magnetic reluctivity [m/H]
+   mu[#{Air,BoundaryNotTerminal}] = muExt;  // permeability [H/m]
+   epsilon[#{Air,BoundaryNotTerminal}] = epsExt;       // permittivity [F/m]
+   sigma[#{Air,BoundaryNotTerminal}]  = sigmaExt;      // conductivity [S/m]
+   
+   
+   // This is a problem with 2 terminals, out of which one is grounded.
+   // The non-grounded can be excited in voltage or in current.
+   // Depending on the excitation type, one of the following functions will be useful.
+   //
+   // When imposing voltage or current via a constraint, the region appears there
+   // If you have a complex, the value has to be a function => use square brackets []
+   Printf(" ************************ Flag_AnalysisTypeFormulation = %g",Flag_AnalysisTypeFormulation);
+   Printf(" ************************ Flag_AnalysisType = %g",Flag_AnalysisType);
+   If (Flag_AnalysisTypeFormulation==0) // E
+	   If((Flag_AnalysisType==0)||(Flag_AnalysisType==1))
+			VTerminal1[] = -VTerm1RMS*Complex[Cos[VTerm1Phase],Sin[VTerm1Phase]];
+	   Else
+			ITerminal1[] = -ITerm1RMS*Complex[Cos[ITerm1Phase],Sin[ITerm1Phase]];
+	   EndIf
+	   If((Flag_AnalysisType==0)||(Flag_AnalysisType==2))
+			VTerminal2[] = -VTerm2RMS*Complex[Cos[VTerm2Phase],Sin[VTerm2Phase]];
+	   Else
+			ITerminal2[] = -ITerm2RMS*Complex[Cos[ITerm2Phase],Sin[ITerm2Phase]];
+	   EndIf
+	Else // H
+	   If((Flag_AnalysisType==0)||(Flag_AnalysisType==1))
+			VTerminal1[] = VTerm1RMS*Complex[Cos[VTerm1Phase],Sin[VTerm1Phase]];
+	   Else
+			ITerminal1[] = ITerm1RMS*Complex[Cos[ITerm1Phase],Sin[ITerm1Phase]];
+	   EndIf
+	   If((Flag_AnalysisType==0)||(Flag_AnalysisType==2))
+			VTerminal2[] = VTerm2RMS*Complex[Cos[VTerm2Phase],Sin[VTerm2Phase]];
+	   Else
+			ITerminal2[] = ITerm2RMS*Complex[Cos[ITerm2Phase],Sin[ITerm2Phase]];
+	   EndIf
+	EndIf
+   
+ }
+
+
+ Constraint{
+
+   // ece BC
+   { Name SetTerminalPotential; Type Assign;  // voltage excited terminals
+     Case {
+       { Region Ground; Value 0.; }
+       If((Flag_AnalysisType==0)||(Flag_AnalysisType==1))
+         { Region Terminal1; Value VTerminal1[]; }
+       EndIf
+	   If((Flag_AnalysisType==0)||(Flag_AnalysisType==2))
+         { Region Terminal2; Value VTerminal2[]; }
+       EndIf
+     }
+   }
+   { Name SetTerminalCurrent; Type Assign;   // current excited terminals
+     Case {
+       If((Flag_AnalysisType==2)|| (Flag_AnalysisType==3))
+         { Region Terminal1; Value ITerminal1[]/h2Ddepth; }  // here the depth is needed
+       EndIf
+	   If((Flag_AnalysisType==1)|| (Flag_AnalysisType==3))
+         { Region Terminal2; Value ITerminal2[]/h2Ddepth; }  // here the depth is needed
+       EndIf
+     }
+   }
+   
+   { Name ImposeE ;
+     Case {
+     }
+   }
+   
+   
+      { Name SetTerminalCurrentH; Type Assign;   // current excited terminals
+     Case {
+       If((Flag_AnalysisType==2)|| (Flag_AnalysisType==3))
+         { Region Cut1Hterminal1; Value ITerminal1[]/h2Ddepth; }  // here the depth is needed
+       EndIf
+	   If((Flag_AnalysisType==1)|| (Flag_AnalysisType==3))
+         { Region Cut1Hterminal2; Value ITerminal2[]/h2Ddepth; }  // here the depth is needed
+       EndIf
+     }
+   }
+   
+   
+  
+  { Name SetTerminalPotentialH; Type Assign;  // voltage excited terminals
+     Case {
+       If((Flag_AnalysisType==0)||(Flag_AnalysisType==1))
+         { Region Cut1Hterminal1; Value VTerminal1[]; }
+       EndIf
+	   If((Flag_AnalysisType==0)||(Flag_AnalysisType==2))
+         { Region Cut1Hterminal2; Value VTerminal2[]; }
+       EndIf
+     }
+   }
+  
+  
+   { Name ImposeHonAxis ;
+     Case {
+       If(Flag_AnalysisTypeFormulation==1) // H
+   	    If(Flag_Axi == 1) // axi
+            { Region Axis; Value 0 ; }
+   		EndIf
+        EndIf
+    }
+   }
+  
+}
+
+modelDir = "../";
+Dir = "res/";
+
+// The formulation and its tools
+
+ If (Flag_AnalysisTypeFormulation==0)  // E inside
+    Include "../formulations/FullWave_E_ece_MIMO2terminals_secondOrder.pro"
+ Else // H inside
+    Include "../formulations/FullWave_H_ece_MIMO2terminals_cplx_2D_and_AXI_secondOrder.pro"
+ EndIf
+ 
+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
+
+PostOperation {
+  If (Flag_AnalysisTypeFormulation==0)   // E
+
+    { Name Maps; NameOfPostProcessing FW_E_ece;
+      Operation {
+      Print [ gradV, OnElementsOf ElementsOf[Vol_FW,OnOneSideOf Sur_FW], File StrCat[Dir,"map_gradV.pos" ]];
+      Print [ E, OnElementsOf Vol_FW, File StrCat[Dir,"map_E.pos" ]];
+      Print [ J, OnElementsOf Vol_FW, File StrCat[Dir,"map_J.pos" ]];
+      Print [ H, OnElementsOf Vol_FW, File StrCat[Dir,"map_H.pos" ]];
+      Print [ Vterminal, OnElementsOf ElementsOf[Sur_FW], File StrCat[Dir,"map_Vterminals.pos" ]];
+      Print [ VnotTerminals, OnElementsOf Vol_FW, File StrCat[Dir,"map_VnotTerminals.pos" ]];
+      Call Change_post_options;
+     }
+    }     
+
+    { Name TransferMatrix; NameOfPostProcessing FW_E_ece ;
+      Operation {
+        If(Flag_AnalysisType==3)  // current excitation
+          Print [ Vterminal, OnRegion Terminal1, Format Table , File  > StrCat[Dir,"test_Z1_RI_formE.txt"]] ;
+		      Print [ Vterminal, OnRegion Terminal2, Format Table , File  > StrCat[Dir,"test_Z2_RI_formE.txt"]] ;  
+        ElseIf(Flag_AnalysisType==0)  // voltage excitation
+		      Print [ Iterminal, OnRegion Terminal1, Format Table , File  > StrCat[Dir,"test_Y1_RI_formE.txt"]] ;
+		      Print [ Iterminal, OnRegion Terminal2, Format Table , File  > StrCat[Dir,"test_Y2_RI_formE.txt"]] ;  
+        ElseIf(Flag_AnalysisType==2)  // hybrid excitation (first ec, second ev)
+          Print [ Vterminal, OnRegion Terminal1, Format Table , File  > StrCat[Dir,"test_H1_RI_formE.txt"]] ;
+		      Print [ Iterminal, OnRegion Terminal2, Format Table , File  > StrCat[Dir,"test_H2_RI_formE.txt"]] ;  
+	      Else // reverse hybrid excitation (first ev, second ec)
+          Print [ Iterminal, OnRegion Terminal1, Format Table , File  > StrCat[Dir,"test_G1_RI_formE.txt"]];
+		      Print [ Vterminal, OnRegion Terminal2, Format Table , File  > StrCat[Dir,"test_G2_RI_formE.txt"]] ;  
+	      EndIf
+      } 
+    }
+  
+  Else // H
+
+    { Name Maps; NameOfPostProcessing FW_H_ece ;
+      Operation {
+      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 [ B, OnElementsOf Vol_FW, File StrCat[Dir,"map_B.pos" ]];
+      Print [ H, OnElementsOf Vol_FW, File StrCat[Dir,"map_H.pos" ]];
+      Print [ rmsH, OnElementsOf Vol_FW, File StrCat[Dir,"map_absH.pos" ]];
+      Call Change_post_options;
+      }
+    }
+
+    { Name TransferMatrix; NameOfPostProcessing FW_H_ece ;
+      Operation {
+        If(Flag_AnalysisType==3)  // current excitation
+          Print [ Vterminal, OnRegion Cut1Hterminal1, Format Table , File  > StrCat[Dir,"test_Z1_RI_formH.txt"]] ;
+          Print [ Vterminal, OnRegion Cut1Hterminal2, Format Table , File  > StrCat[Dir,"test_Z2_RI_formH.txt"]] ;   
+        ElseIf(Flag_AnalysisType==0)  // voltage excitation
+          Print [ Iterminal, OnRegion Cut1Hterminal1, Format Table , File  > StrCat[Dir,"test_Y1_RI_formH.txt"]] ;
+          Print [ Iterminal, OnRegion Cut1Hterminal2, Format Table , File  > StrCat[Dir,"test_Y2_RI_formH.txt"]] ;  
+        ElseIf(Flag_AnalysisType==2)  // hybrid excitation (first ec, second ev)
+          Print [ Vterminal, OnRegion Cut1Hterminal1, Format Table , File  > StrCat[Dir,"test_H1_RI_formH.txt"]] ;
+          Print [ Iterminal, OnRegion Cut1Hterminal2, Format Table , File  > StrCat[Dir,"test_H2_RI_formH.txt"]] ;  
+        Else // reverse hybrid excitation (first ev, second ec)
+          Print [ Iterminal, OnRegion Cut1Hterminal1, Format Table , File  > StrCat[Dir,"test_G1_RI_formH.txt"]] ;
+          Print [ Vterminal, OnRegion Cut1Hterminal2, Format Table , File  > StrCat[Dir,"test_G2_RI_formH.txt"]] ;  
+        EndIf
+      }  
+    }
+ 
+  EndIf
+
+}
diff --git a/cylinderCoax/geo_pro/coax2Daxi_data.pro b/cylinderCoax/geo_pro/coax2Daxi_data.pro
new file mode 100644
index 0000000..a7d8854
--- /dev/null
+++ b/cylinderCoax/geo_pro/coax2Daxi_data.pro
@@ -0,0 +1,178 @@
+colorro = "LightGrey";
+
+// new menu for the interface - holds geometrical data
+cGUI= "{TEST-Coax2Daxi-ProblemData/";
+mGeo = StrCat[cGUI,"0Geometrical parameters/0"];
+colorMGeo = "Ivory";
+close_menu = 1;
+DefineConstant[
+  a = {4e-3, Name StrCat[mGeo, "0a=radius of inner conductor"], Units "m", Highlight Str[colorMGeo], Closed close_menu},
+  b = {8e-3, Name StrCat[mGeo, "0a=radius of outer conductor"], Units "m", Highlight Str[colorMGeo], Closed close_menu},
+  l = {1, Name StrCat[mGeo, "1l=cable length"], Units "m", Highlight Str[colorMGeo]}
+];
+
+// new menu for the interface - holds material data
+mMat = StrCat[cGUI,"1Material parameters/0"];
+colorMMat = "AliceBlue";
+DefineConstant[
+  epsrInt = {1, Name StrCat[mMat, "0Relative permittivity int part"], Units "-", Highlight Str[colorMMat], Closed close_menu},
+  murInt = {1, Name StrCat[mMat, "1Relative permeability int part"], Units "-", Highlight Str[colorMMat]},
+  sigmaInt = {6.6e7, Name StrCat[mMat, "2Conductivity int part"], Units "S/m", Highlight Str[colorMMat]},
+  epsrExt = {1, Name StrCat[mMat, "3Relative permittivity ext part"], Units "-", Highlight Str[colorMMat], Closed close_menu},
+  murExt = {1, Name StrCat[mMat, "4Relative permeability ext part"], Units "-", Highlight Str[colorMMat]},
+  sigmaExt = {0, Name StrCat[mMat, "5Conductivity ext part"], Units "S/m", Highlight Str[colorMMat]}
+];
+
+eps0 = 8.854187818e-12; // F/m - absolute permitivity of vacuum
+mu0  = 4.e-7 * Pi;      // H/m - absolute permeability of vacuum
+nu0  = 1/mu0;           // m/H - absolute reluctivity of vacuum
+
+epsInt = eps0*epsrInt;        // F/m - absolute permitivity  - int
+muInt = mu0*murInt;           // F/m - absolute permeability - int
+nuInt = 1/muInt;              // m/H - absolute reluctivity  - int
+
+epsExt = eps0*epsrExt;        // F/m - absolute permitivity  - ext
+muExt = mu0*murExt;           // F/m - absolute permeability - ext
+nuExt = 1/muExt;              // m/H - absolute reluctivity  - ext
+
+// new menu for the interface - type of BC
+mTypeBC = StrCat[cGUI,"2Type of BC/0"];
+colorMTypeBC = "Red";
+
+// new menu for the interface - type of formulation
+mTypeBC2 = StrCat[cGUI,"3Type of formulation/0"];
+colorMTypeBC = "Blue";
+DefineConstant[
+  Flag_AnalysisTypeFormulation = { FORMULATION,
+    Choices{
+      0="0: formulation in E inside + ECE",
+      1="1: formulation in H inside + ECE"
+    },
+    Name Str[mTypeBC2], Highlight Str[colorMTypeBC], Closed !close_menu,
+    Help Str["E inside, electric V on the boundary",
+      "H inside, reduced magnetic V on the boundary"],
+    ServerAction Str["Reset", "GetDP/1ResolutionChoices"]}
+];
+
+// new menu for the interface - values of excitations - significance depends on the type of BC
+mValuesBC = StrCat[cGUI,"4Values of BC (frequency analysis)/0"];
+colorMValuesBC = "Ivory";
+
+If (MODEL_ID2 == 1) // HF simulations
+	fmin = 1e7; fmax = 6e8; nop = 100; 
+	freqs()=LinSpace[fmin,fmax,nop];
+EndIf
+If  (MODEL_ID2 == 2) // LF simulations
+	fmin = 1;   fmax = 1e7; nop = 20;
+	freqs()=LogSpace[Log10[fmin],Log10[fmax],nop];
+EndIf
+
+DefineConstant[
+  _use_freq_loop = {0,  Choices{0,1}, Name StrCat[mValuesBC, "0Loop on frequencies?"],
+    ServerAction Str["Reset", StrCat[mValuesBC, "0Working Frequency"]], Highlight "Yellow"}
+  Freq  = {freqs(0), Choices{freqs()}, Loop _use_freq_loop, Name StrCat[mValuesBC, "0Working Frequency"],
+    Units "Hz", Highlight  Str[colorMValuesBC], Closed !close_menu }
+];
+
+NbTerminals = 2; // No of terminals which are not grounded
+
+Type1 = (TERMINAL_EXC == 0)  ? 2 : 1 ;  // 2 ==> voltage exc
+Type2 = (TERMINAL_EXC2 == 0) ? 2 : 1 ; 
+	
+DefineConstant[
+	TypeExcTerminal1 = {Type1,Choices{1="ec",2="ev"},Name StrCat[mValuesBC, "4Exitation type/0Terminal 1 (bottom)"]}
+	TypeExcTerminal2 = {Type2,Choices{1="ec",2="ev"},Name StrCat[mValuesBC, "4Exitation type/1Terminal 2 (top)"]}
+	ActiveTerminal = {1,Choices{1="1",2="2"},Name StrCat[mValuesBC, "4Exitation type/Active terminal"]}
+
+	VGroundRMS  = {0, Name StrCat[mValuesBC, "3Ground-CylSurf/1Potential on gnd (rms)"], Units "V" , 
+		ReadOnly 1, Highlight Str[colorro]}
+  	VGroundPhase  = {0, Name StrCat[mValuesBC, "3Ground-CylSurf/1Potential on gnd, phase"], Units "rad",  
+		ReadOnly 1, Highlight Str[colorro]}
+];
+
+If (ActiveTerminal == 1)
+	PassiveTerminal = 2;
+Else
+	PassiveTerminal = 1;
+EndIf
+  
+If ((TypeExcTerminal1==1)&(TypeExcTerminal2==1))  // current exc
+	Flag_AnalysisType = 3;
+	Printf("======================> Both terminals are current excited => Z transfer matrix");
+	Printf("======================> Active terminal is %g",ActiveTerminal);
+	Printf("======================> Two s1p files are created, containing Z%g%g and Z%g%g",ActiveTerminal,ActiveTerminal,PassiveTerminal,ActiveTerminal);
+	If (ActiveTerminal == 1)
+		ITerm1RMS = 1; ITerm1Phase = 0;
+		ITerm2RMS = 0; ITerm2Phase = 0;
+		If (Flag_AnalysisTypeFormulation==1)   // H - to understand why I have to do this
+		   ITerm1RMS = -1;
+		EndIf
+	Else
+	    ITerm1RMS = 0; ITerm1Phase = 0;
+		ITerm2RMS = 1; ITerm2Phase = 0;
+	EndIf
+ElseIf ((TypeExcTerminal1==2)&(TypeExcTerminal2==2))  // voltage exc
+	Flag_AnalysisType = 0;
+	Printf("======================> Both terminals are voltage excited => Y transfer matrix");
+	Printf("======================> Active terminal is %g",ActiveTerminal);
+	Printf("======================> Two s1p files are created, containing Y%g%g and Y%g%g",ActiveTerminal,ActiveTerminal,PassiveTerminal,ActiveTerminal);
+	If (ActiveTerminal == 1)
+		VTerm1RMS = 1; VTerm1Phase = 0;
+		VTerm2RMS = 0; VTerm2Phase = 0;
+	Else
+	    VTerm1RMS = 0; VTerm1Phase = 0;
+		VTerm2RMS = (Flag_AnalysisTypeFormulation==1) ? -1 : 1; 
+		// H - constraint is linked to an automatically generated cut, which is oriented and depends on geo/msh
+		// the change of sign here is due to the orientation of this cut. It may need to be adapted with another geo/msh
+		VTerm2Phase = 0;  
+	EndIf
+ElseIf ((TypeExcTerminal1==1)&(TypeExcTerminal2==2))  //hybrid H exc
+	Flag_AnalysisType = 2;
+	Printf("======================> Terminal 1 (top) is current excited, Terminal 2 (right) is voltage excited => H transfer matrix");
+	Printf("======================> Active terminal is %g",ActiveTerminal);
+	Printf("======================> Two s1p files are created, containing H%g%g and H%g%g",ActiveTerminal,ActiveTerminal,PassiveTerminal,ActiveTerminal);
+	If (ActiveTerminal == 1)
+		ITerm1RMS = (Flag_AnalysisTypeFormulation==1)?-1:1; // change of sign due to cut orientation
+		ITerm1Phase = 0;
+		VTerm2RMS = 0; VTerm2Phase = 0;
+	Else
+	    ITerm1RMS = 0; ITerm1Phase = 0;
+		VTerm2RMS = 1; VTerm2Phase = 0;  
+	EndIf
+Else  //hybrid G exc
+	Flag_AnalysisType = 1;
+	Printf("======================> Terminal 1 (top) is voltage excited, Terminal 2 (right) is current excited => G transfer matrix");
+	Printf("======================> Active terminal is %g",ActiveTerminal);
+	Printf("======================> Two s1p files are created, containing G%g%g and G%g%g",ActiveTerminal,ActiveTerminal,PassiveTerminal,ActiveTerminal);
+	If (ActiveTerminal == 1)
+		VTerm1RMS = 1; VTerm1Phase = 0;
+		ITerm2RMS = 0; ITerm2Phase = 0;
+	Else
+	    VTerm1RMS = 0; VTerm1Phase = 0;
+		ITerm2RMS = 1; ITerm2Phase = 0;
+	EndIf
+EndIf
+
+meshSettings = StrCat[cGUI,"5MeshSettings/0"];
+
+DefineConstant[
+  _use_transfinite = {1, Choices{0,1}, Name StrCat[meshSettings, "0Transfinite mesh?"], Highlight "Yellow", Closed !close_menu}
+   s = {1, Name  StrCat[meshSettings, "1Mesh factor"]}
+   Flag_ElementOrder = { FE_ORDER,
+    Choices{
+      1="1: first order elements",
+      2="2: second order elements"
+    },
+	Name StrCat[meshSettings, "2Element order"], Highlight "Red", Closed !close_menu,
+    ServerAction Str["Reset", "GetDP/1ResolutionChoices"]}
+	_use_recombine = {0, Choices{0,1}, Name StrCat[meshSettings, "3Use recombine?"], Highlight "Yellow", Closed !close_menu}
+];
+
+FE_Order = Flag_ElementOrder;
+
+Printf("FE_Order = %g",FE_Order);
+
+modelDim = 2; 
+Flag_Axi = 1; 
+h2Ddepth = (modelDim==3) ? 1 : (Flag_Axi ? 2*Pi : l);
+
diff --git a/cylinderCoax/geo_pro/coax3D.geo b/cylinderCoax/geo_pro/coax3D.geo
new file mode 100644
index 0000000..cddc159
--- /dev/null
+++ b/cylinderCoax/geo_pro/coax3D.geo
@@ -0,0 +1,99 @@
+Include "coax3D_data.pro";
+
+p =0.25; // portion of a
+
+nop = 4/4;  // no of slices per 1/4 circle
+nopR = s*4; //5; //5*2; // no of points on the Oa and atob points
+nopA = s*5; //5*4;
+
+
+pc = newp; Point(pc) = {0,0,0};
+angle = Pi/2/nop;
+
+For k1 In {0:4*nop-1}
+  alpha = angle*k1;
+  pa[] += newp ; Point(pa[k1]) = {a*Cos[alpha],a*Sin[alpha],0};
+
+  // for controlling the mesh a bit better
+  pp[] += newp ; Point(pp[k1]) = {p*a*Cos[alpha],p*a*Sin[alpha],0};
+  la[] += newl ; Line(la[k1]) = {pp[k1],pa[k1]};
+
+  pb[] += newp ; Point(pb[k1]) = {b*Cos[alpha],b*Sin[alpha],0};
+  lb[] += newl ; Line(lb[k1]) = {pa[k1],pb[k1]};
+EndFor
+
+For k1 In {0:4*nop-1}
+  k2 = (k1==4*nop-1) ? 0 : k1+1;
+  cla[] += newl ; Circle(newl) = {pa[k1],pc,pa[k2]};
+  clp[] += newl ; Circle(newl) = {pp[k1],pc,pp[k2]}; 
+  clb[] += newl ; Circle(newl) = {pb[k1],pc,pb[k2]};
+
+  //In
+  Curve Loop(newcl) = {la[k1],cla[k1],-la[k2], -clp[k1]};
+  scable_a[] += news; Plane Surface(news) = newcl-1;
+  // Out
+  Curve Loop(newcl) = {lb[k1],clb[k1],-lb[k2],-cla[k1]};
+  scable_b[] += news; Plane Surface(news) = newcl-1;
+EndFor
+Curve Loop(newcl) = {clp[]};
+scable_a[] += news; Plane Surface(news) = newcl-1; // center
+
+Transfinite Curve{la[]} = nopR; 
+Transfinite Curve{lb[]} = nopR;
+
+Transfinite Curve{cla[],clb[],clp[]} = nopA;
+Transfinite Surface{scable_a[]};
+Transfinite Surface{scable_b[]};
+
+If (_use_recombine)
+  Recombine Surface{scable_a[]};
+  Recombine Surface{scable_b[]};
+EndIf
+
+nbPointAlongOz = Ceil[l/lcarLambda];
+Printf("nbPointAlongOz = %g",nbPointAlongOz);
+//nbPointAlongOz = 16;  // just to see if I can generate second order systems
+
+// You may use a loop, much easier
+For k1 In {0:4*nop}
+  aux_a[]=Extrude {0,0,l} {
+    Surface{scable_a[k1]}; Layers{ {nbPointAlongOz}, {1} }; Recombine _use_recombine;
+  }; // array containing: 0=top surface; 1=volume; the rest of the surfaces
+  vcable_a[] += aux_a[1]; // volume
+  scable_aa[] += aux_a[0]; // top surface after extruding
+EndFor
+
+For k1 In {0:4*nop-1}
+  aux_b[]=Extrude {0,0,l} {
+    Surface{scable_b[k1]}; Layers{ {nbPointAlongOz}, {1} }; Recombine _use_recombine;
+  };
+  vcable_b[] += aux_b[1]; // volume
+  scable_bb[] += aux_b[0]; // top surface after extruding
+EndFor
+
+// and no need of hunting :-)
+ground[] = Abs[CombinedBoundary{Volume{:};}] ; // Total boundary + getting rid of possible 'signed' surfaces
+ground[] -= {scable_a[],scable_aa[], scable_b[], scable_bb[]}; // Removing top and bottom from array
+
+
+
+// Physical regions
+// --------------------------
+
+Physical Volume ("Conductor", 100) = {vcable_a[]} ;         
+Physical Volume ("Air", 101) = {vcable_b[]} ;
+
+Physical Surface ("Ground",   120) = {ground[]} ;          
+Physical Surface ("Terminal1", 121) = {scable_a[]} ;       
+Physical Surface ("Terminal2", 122) = {scable_aa[]} ;      
+Physical Surface ("BottomBoundary", 133) = {scable_b[]} ;
+Physical Surface ("TopBoundary", 134) = {scable_bb[]} ;
+
+Physical Surface ("BoundaryNotTerminal", 131) = {scable_b[],scable_bb[]} ;
+
+// cut for the H formulation
+Cohomology(1) {{131},{}}; 
+
+// For aestetics
+Recursive Color SkyBlue { Volume{vcable_b[]};}
+Recursive Color Red  { Volume{vcable_a[]};}
diff --git a/cylinderCoax/geo_pro/coax3D.pro b/cylinderCoax/geo_pro/coax3D.pro
new file mode 100644
index 0000000..7635b39
--- /dev/null
+++ b/cylinderCoax/geo_pro/coax3D.pro
@@ -0,0 +1,217 @@
+Include "coax3D_data.pro";
+
+Group{
+  Conductor = Region[100]; 
+  Air = Region[101];
+
+  Terminal1 = Region[121]; 
+  Terminal2 = Region[122]; 
+
+  Ground = Region[120];    
+  BottomBoundary = Region[133];
+  TopBoundary = Region[134];
+
+
+  BoundaryNotTerminal = Region[{BottomBoundary,TopBoundary}];
+  Sur_Terminals_FWece = Region[{Terminal1,Terminal2,Ground}];
+
+  Vol_FW =  Region[{Conductor,Air}];
+
+  Sur_FW = Region[{BoundaryNotTerminal,Sur_Terminals_FWece}]; //all boundary
+  Dom_FW = Region[ {Vol_FW, Sur_FW} ];
+
+  Cut1Hterminal1 = Region[{135}];
+  Cut1Hterminal2 = Region[{136}];
+  Cut1H = Region[{Cut1Hterminal1,Cut1Hterminal2}];
+}
+
+Function{
+  // Material properties are defined piecewise, for elementary groups
+  // the values in the right hand side are defined in Ishape2d_data.pro
+  
+  nu[#{Conductor,Terminal1,Terminal2}] = nuInt;  // magnetic reluctivity [m/H]
+  mu[#{Conductor,Terminal1,Terminal2}] = muInt;  // permeability [H/m]
+  epsilon[#{Conductor,Terminal1,Terminal2}] = epsInt;       // permittivity [F/m]
+  sigma[#{Conductor,Terminal1,Terminal2}]  = sigmaInt;      // conductivity [S/m]
+
+  nu[#{Air,BoundaryNotTerminal,Ground}] = nuExt;  // magnetic reluctivity [m/H]
+  mu[#{Air,BoundaryNotTerminal,Ground}] = muExt;  // permeability [H/m]
+  epsilon[#{Air,BoundaryNotTerminal,Ground}] = epsExt;       // permittivity [F/m]
+  sigma[#{Air,BoundaryNotTerminal,Ground}]  = sigmaExt;      // conductivity [S/m]
+  
+  // This is a problem with 2 terminals, out of which one is grounded.
+  // The non-grounded can be excited in voltage or in current.
+  // Depending on the excitation type, one of the following functions will be useful.
+  //
+  // When imposing voltage or current via a constraint, the region appears there
+  // If you have a complex, the value has to be a function => use square brackets []
+  If (Flag_AnalysisTypeFormulation==0) // E
+    If((Flag_AnalysisType==0)||(Flag_AnalysisType==1))
+    VTerminal1[] = -VTerm1RMS*Complex[Cos[VTerm1Phase],Sin[VTerm1Phase]];
+    Else
+    ITerminal1[] = -ITerm1RMS*Complex[Cos[ITerm1Phase],Sin[ITerm1Phase]];
+    EndIf
+    If((Flag_AnalysisType==0)||(Flag_AnalysisType==2))
+    VTerminal2[] = -VTerm2RMS*Complex[Cos[VTerm2Phase],Sin[VTerm2Phase]];
+    Else
+    ITerminal2[] = -ITerm2RMS*Complex[Cos[ITerm2Phase],Sin[ITerm2Phase]];
+    EndIf
+	Else // H
+    If((Flag_AnalysisType==0)||(Flag_AnalysisType==1))
+    VTerminal1[] = VTerm1RMS*Complex[Cos[VTerm1Phase],Sin[VTerm1Phase]];
+    Else
+    ITerminal1[] = ITerm1RMS*Complex[Cos[ITerm1Phase],Sin[ITerm1Phase]];
+    EndIf
+    If((Flag_AnalysisType==0)||(Flag_AnalysisType==2))
+    VTerminal2[] = VTerm2RMS*Complex[Cos[VTerm2Phase],Sin[VTerm2Phase]];
+    Else
+    ITerminal2[] = ITerm2RMS*Complex[Cos[ITerm2Phase],Sin[ITerm2Phase]];
+    EndIf
+	EndIf  
+}
+
+Constraint{
+  // ece BC
+  { Name SetTerminalPotential; Type Assign;  // voltage excited terminals
+    Case {
+      { Region Ground; Value 0.; }
+      If((Flag_AnalysisType==0)||(Flag_AnalysisType==1))
+        { Region Terminal1; Value VTerminal1[]; }
+      EndIf
+      If((Flag_AnalysisType==0)||(Flag_AnalysisType==2))
+        { Region Terminal2; Value VTerminal2[]; }
+      EndIf
+    }
+  }
+   
+  { Name SetTerminalCurrent; Type Assign;   // current excited terminals
+    Case {
+      If((Flag_AnalysisType==2)|| (Flag_AnalysisType==3))
+        { Region Terminal1; Value ITerminal1[]/h2Ddepth; }  // depth is needed
+      EndIf
+      If((Flag_AnalysisType==1)|| (Flag_AnalysisType==3))
+        { Region Terminal2; Value ITerminal2[]/h2Ddepth; }  // depth is needed
+      EndIf
+    }
+   }
+   
+  { Name ImposeE ;
+    Case {
+    }
+  }
+  
+  { Name SetTerminalCurrentH; Type Assign;   // current excited terminals
+    Case {
+      If((Flag_AnalysisType==2)|| (Flag_AnalysisType==3))
+      { Region Cut1Hterminal1; Value ITerminal1[]/h2Ddepth; }  // here the depth is needed
+      EndIf
+	    If((Flag_AnalysisType==1)|| (Flag_AnalysisType==3))
+        { Region Cut1Hterminal2; Value ITerminal2[]/h2Ddepth; }  // here the depth is needed
+      EndIf
+    }
+  }
+   
+  { Name SetTerminalPotentialH; Type Assign;  // voltage excited terminals
+    Case {
+      If((Flag_AnalysisType==0)||(Flag_AnalysisType==1))
+        { Region Cut1Hterminal1; Value VTerminal1[]; }
+      EndIf
+      If((Flag_AnalysisType==0)||(Flag_AnalysisType==2))
+        { Region Cut1Hterminal2; Value VTerminal2[]; }
+      EndIf
+    }
+  }
+  
+}
+
+modelDir = "../";
+Dir = "res/";
+ 
+// The formulation and its tools
+If (Flag_AnalysisTypeFormulation==0)  // E inside
+  Include "../formulations/FullWave_E_ece_MIMO2terminals_secondOrder.pro"
+Else // H inside
+  Include "../formulations/FullWave_H_ece_MIMO2terminals_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
+
+PostOperation {
+
+  If (Flag_AnalysisTypeFormulation==0)   // E
+
+    { Name Maps; NameOfPostProcessing FW_E_ece;
+      Operation {
+      Print [ gradV, OnElementsOf ElementsOf[Vol_FW,OnOneSideOf Sur_FW], File StrCat[Dir,"map_gradV.pos" ]];
+      Print [ E, OnElementsOf Vol_FW, File StrCat[Dir,"map_E.pos" ]];
+      Print [ J, OnElementsOf Vol_FW, File StrCat[Dir,"map_J.pos" ]];
+      Print [ H, OnElementsOf Vol_FW, File StrCat[Dir,"map_H.pos" ]];
+      Print [ Vterminal, OnElementsOf ElementsOf[Sur_FW], File StrCat[Dir,"map_Vterminals.pos" ]];
+      Print [ VnotTerminals, OnElementsOf Vol_FW, File StrCat[Dir,"map_VnotTerminals.pos" ]];
+      Call Change_post_options;
+     }
+    }     
+
+	  { Name TransferMatrix; NameOfPostProcessing FW_E_ece ;
+      Operation {
+        If(Flag_AnalysisType==3)  // current excitation
+          Print [ Vterminal, OnRegion Terminal1, Format Table , File  > StrCat[Dir,"test_Z1_RI_formE.txt"]] ;
+		      Print [ Vterminal, OnRegion Terminal2, Format Table , File  > StrCat[Dir,"test_Z2_RI_formE.txt"]] ;  
+        ElseIf(Flag_AnalysisType==0)  // voltage excitation
+		      Print [ Iterminal, OnRegion Terminal1, Format Table , File  > StrCat[Dir,"test_Y1_RI_formE.txt"]] ;
+		      Print [ Iterminal, OnRegion Terminal2, Format Table , File  > StrCat[Dir,"test_Y2_RI_formE.txt"]] ;  
+        ElseIf(Flag_AnalysisType==2)  // hybrid excitation (first ec, second ev)
+          Print [ Vterminal, OnRegion Terminal1, Format Table , File  > StrCat[Dir,"test_H1_RI_formE.txt"]] ;
+		      Print [ Iterminal, OnRegion Terminal2, Format Table , File  > StrCat[Dir,"test_H2_RI_formE.txt"]] ;  
+	      Else // reverse hybrid excitation (first ev, second ec)
+          Print [ Iterminal, OnRegion Terminal1, Format Table , File  > StrCat[Dir,"test_G1_RI_formE.txt"]] ;
+		      Print [ Vterminal, OnRegion Terminal2, Format Table , File  > StrCat[Dir,"test_G2_RI_formE.txt"]] ;  
+	      EndIf
+      } 
+    }
+	
+  Else // H
+
+    { Name Maps; NameOfPostProcessing FW_H_ece ;
+      Operation {
+      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 [ B, OnElementsOf Vol_FW, File StrCat[Dir,"map_B.pos" ]];
+      Print [ H, OnElementsOf Vol_FW, File StrCat[Dir,"map_H.pos" ]];
+      Print [ rmsH, OnElementsOf Vol_FW, File StrCat[Dir,"map_absH.pos" ]];
+      Call Change_post_options;
+      }
+    }
+
+    { Name TransferMatrix; NameOfPostProcessing FW_H_ece ;
+      Operation {
+        If(Flag_AnalysisType==3)  // current excitation
+          Print [ Vterminal, OnRegion Cut1Hterminal1, Format Table , File  > StrCat[Dir,"test_Z1_RI_formH.txt"]] ;
+          Print [ Vterminal, OnRegion Cut1Hterminal2, Format Table , File  > StrCat[Dir,"test_Z2_RI_formH.txt"]] ;   
+        ElseIf(Flag_AnalysisType==0)  // voltage excitation
+          Print [ Iterminal, OnRegion Cut1Hterminal1, Format Table , File  > StrCat[Dir,"test_Y1_RI_formH.txt"]] ;
+          Print [ Iterminal, OnRegion Cut1Hterminal2, Format Table , File  > StrCat[Dir,"test_Y2_RI_formH.txt"]] ;  
+        ElseIf(Flag_AnalysisType==2)  // hybrid excitation (first ec, second ev)
+          Print [ Vterminal, OnRegion Cut1Hterminal1, Format Table , File  > StrCat[Dir,"test_H1_RI_formH.txt"]] ;
+          Print [ Iterminal, OnRegion Cut1Hterminal2, Format Table , File  > StrCat[Dir,"test_H2_RI_formH.txt"]] ;  
+        Else // reverse hybrid excitation (first ev, second ec)
+          Print [ Iterminal, OnRegion Cut1Hterminal1, Format Table , File  > StrCat[Dir,"test_G1_RI_formH.txt"]] ;
+          Print [ Vterminal, OnRegion Cut1Hterminal2, Format Table , File  > StrCat[Dir,"test_G2_RI_formH.txt"]] ;  
+        EndIf
+      }  
+    }
+      
+  EndIf
+
+}
diff --git a/cylinderCoax/geo_pro/coax3D_data.pro b/cylinderCoax/geo_pro/coax3D_data.pro
new file mode 100644
index 0000000..6b11b75
--- /dev/null
+++ b/cylinderCoax/geo_pro/coax3D_data.pro
@@ -0,0 +1,174 @@
+
+
+colorro = "LightGrey";
+
+// new menu for the interface - holds geometrical data
+cGUI= "{TEST-Coax3D-ProblemData/";
+mGeo = StrCat[cGUI,"0Geometrical parameters/0"];
+colorMGeo = "Ivory";
+close_menu = 1;
+DefineConstant[
+  a = {4e-3, Name StrCat[mGeo, "0a=radius of inner conductor"], Units "m", Highlight Str[colorMGeo], Closed close_menu},
+  b = {8e-3, Name StrCat[mGeo, "0a=radius of outer conductor"], Units "m", Highlight Str[colorMGeo], Closed close_menu},
+  l = {1, Name StrCat[mGeo, "1l=cable length"], Units "m", Highlight Str[colorMGeo]}
+];
+
+/* new menu for the interface - holds material data */
+mMat = StrCat[cGUI,"1Material parameters/0"];
+colorMMat = "AliceBlue";
+DefineConstant[
+  epsrInt = {1, Name StrCat[mMat, "0Relative permittivity int part"], Units "-", Highlight Str[colorMMat], Closed close_menu},
+  murInt = {1, Name StrCat[mMat, "1Relative permeability int part"], Units "-", Highlight Str[colorMMat]},
+  sigmaInt = {6.6e7, Name StrCat[mMat, "2Conductivity int part"], Units "S/m", Highlight Str[colorMMat]},
+  epsrExt = {1, Name StrCat[mMat, "3Relative permittivity ext part"], Units "-", Highlight Str[colorMMat], Closed close_menu},
+  murExt = {1, Name StrCat[mMat, "4Relative permeability ext part"], Units "-", Highlight Str[colorMMat]},
+  sigmaExt = {0, Name StrCat[mMat, "5Conductivity ext part"], Units "S/m", Highlight Str[colorMMat]}
+];
+
+eps0 = 8.854187818e-12; // F/m - absolute permitivity of vacuum
+mu0  = 4.e-7 * Pi;      // H/m - absolute permeability of vacuum
+nu0  = 1/mu0;           // m/H - absolute reluctivity of vacuum
+
+epsInt = eps0*epsrInt;        // F/m - absolute permitivity  - int
+muInt = mu0*murInt;           // F/m - absolute permeability - int
+nuInt = 1/muInt;              // m/H - absolute reluctivity  - int
+//
+epsExt = eps0*epsrExt;        // F/m - absolute permitivity  - ext
+muExt = mu0*murExt;           // F/m - absolute permeability - ext
+nuExt = 1/muExt;              // m/H - absolute reluctivity  - ext
+
+// new menu for the interface - type of BC 
+mTypeBC = StrCat[cGUI,"2Type of BC/0"];
+colorMTypeBC = "Red";
+
+
+// new menu for the interface - type of formulation
+mTypeBC2 = StrCat[cGUI,"3Type of formulation/0"];
+colorMTypeBC = "Blue";
+DefineConstant[
+  Flag_AnalysisTypeFormulation = { 0,
+    Choices{
+      0="0: formulation in E inside + ECE",
+      1="1: formulation in H inside + ECE"
+    },
+    Name Str[mTypeBC2], Highlight Str[colorMTypeBC], Closed !close_menu,
+    Help Str["E inside, electric V on the boundary",
+      "H inside, reduced magnetic V on the boundary"],
+    ServerAction Str["Reset", "GetDP/1ResolutionChoices"]}
+];
+
+// new menu for the interface - values of excitations - significance depends on the type of BC
+mValuesBC = StrCat[cGUI,"4Values of BC (frequency analysis)/0"];
+colorMValuesBC = "Ivory";
+
+fmin = 1e7;
+fmin = 3e8;      
+fmax = 6e8; 
+nop = 50; 
+
+freqs()=LinSpace[fmin,fmax,nop];
+
+DefineConstant[
+  _use_freq_loop = {0,  Choices{0,1}, Name StrCat[mValuesBC, "0Loop on frequencies?"],
+    ServerAction Str["Reset", StrCat[mValuesBC, "0Working Frequency"]], Highlight "Yellow"}
+    Freq  = {freqs(0), Choices{freqs()}, Loop _use_freq_loop, Name StrCat[mValuesBC, "0Working Frequency"],
+    Units "Hz", Highlight  Str[colorMValuesBC], Closed !close_menu }
+];
+
+
+NbTerminals = 2; // No of terminals which are not grounded
+
+DefineConstant[
+  TypeExcTerminal1 = {2,Choices{1="ec",2="ev"},Name StrCat[mValuesBC, "4Exitation type/0Terminal 1 (bottom)"]}
+  TypeExcTerminal2 = {2,Choices{1="ec",2="ev"},Name StrCat[mValuesBC, "4Exitation type/1Terminal 2 (top)"]}
+  ActiveTerminal = {1,Choices{1="1",2="2"},Name StrCat[mValuesBC, "4Exitation type/Active terminal"]}
+
+  VGroundRMS  = {0, Name StrCat[mValuesBC, "3Ground-CylSurf/1Potential on gnd (rms)"], Units "V" , 
+    ReadOnly 1, Highlight Str[colorro]}
+  VGroundPhase  = {0, Name StrCat[mValuesBC, "3Ground-CylSurf/1Potential on gnd, phase"], Units "rad",  
+    ReadOnly 1, Highlight Str[colorro]}
+];
+
+
+PassiveTerminal = (ActiveTerminal == 1)?2:1;
+
+meshSettings = StrCat[cGUI,"5MeshSettings/0"];
+
+DefineConstant[
+  s = {1, Name StrCat[meshSettings, "1Mesh factor"]}
+  Flag_ElementOrder = { 1,
+  Choices{
+    1="1: first order elements",
+    2="2: second order elements"
+  },
+  Name StrCat[meshSettings, "2Element order"], Highlight "Red", Closed !close_menu,
+  ServerAction Str["Reset", "GetDP/1ResolutionChoices"]}
+	_use_recombine = {0, Choices{0,1}, Name StrCat[meshSettings, "3Use recombine?"], Highlight "Yellow", Closed !close_menu}
+];
+
+FE_Order = Flag_ElementOrder;
+
+c0 = 3e8;
+lambdaMin = c0/fmax;
+lcarLambda = s*lambdaMin/10;
+
+If ((TypeExcTerminal1==1)&(TypeExcTerminal2==1))  // current exc
+	Flag_AnalysisType = 3;
+	Printf("======================> Both terminals are current excited => Z transfer matrix");
+	Printf("======================> Active terminal is %g",ActiveTerminal);
+	Printf("======================> Two s1p files are created, containing Z%g%g and Z%g%g",ActiveTerminal,ActiveTerminal,PassiveTerminal,ActiveTerminal);
+	If (ActiveTerminal == 1)
+		ITerm1RMS = 1; ITerm1Phase = 0;
+		ITerm2RMS = 0; ITerm2Phase = 0;
+	Else
+	    ITerm1RMS = 0; ITerm1Phase = 0;
+		ITerm2RMS = 1; ITerm2Phase = 0;
+	EndIf
+ElseIf ((TypeExcTerminal1==2)&(TypeExcTerminal2==2))  // voltage exc
+	Flag_AnalysisType = 0;
+	Printf("======================> Both terminals are voltage excited => Y transfer matrix");
+	Printf("======================> Active terminal is %g",ActiveTerminal);
+	Printf("======================> Two s1p files are created, containing Y%g%g and Y%g%g",ActiveTerminal,ActiveTerminal,PassiveTerminal,ActiveTerminal);
+	If (ActiveTerminal == 1)
+		VTerm1RMS = 1; VTerm1Phase = 0;
+		VTerm2RMS = 0; VTerm2Phase = 0;
+	Else
+	    VTerm1RMS = 0; VTerm1Phase = 0;
+		VTerm2RMS = 1; VTerm2Phase = 0;  
+		If(Flag_AnalysisTypeFormulation==1)  // H - uncomment when no recombine
+		  If(_use_recombine==0)  
+		   VTerm2RMS = -1;
+		  EndIf
+		EndIf
+	EndIf
+ElseIf ((TypeExcTerminal1==1)&(TypeExcTerminal2==2))  //hybrid H exc
+	Flag_AnalysisType = 2;
+	Printf("======================> Terminal 1 (top) is current excited, Terminal 2 (right) is voltage excited => H transfer matrix");
+	Printf("======================> Active terminal is %g",ActiveTerminal);
+	Printf("======================> Two s1p files are created, containing H%g%g and H%g%g",ActiveTerminal,ActiveTerminal,PassiveTerminal,ActiveTerminal);
+	If (ActiveTerminal == 1)
+		ITerm1RMS = 1; ITerm1Phase = 0;
+		VTerm2RMS = 0; VTerm2Phase = 0;
+	Else
+	    ITerm1RMS = 0; ITerm1Phase = 0;
+		VTerm2RMS = 1; VTerm2Phase = 0;  
+	EndIf
+Else  //hybrid G exc
+	Flag_AnalysisType = 1;
+	Printf("======================> Terminal 1 (top) is voltage excited, Terminal 2 (right) is current excited => G transfer matrix");
+	Printf("======================> Active terminal is %g",ActiveTerminal);
+	Printf("======================> Two s1p files are created, containing G%g%g and G%g%g",ActiveTerminal,ActiveTerminal,PassiveTerminal,ActiveTerminal);
+	If (ActiveTerminal == 1)
+		VTerm1RMS = 1; VTerm1Phase = 0;
+		ITerm2RMS = 0; ITerm2Phase = 0;
+	Else
+	    VTerm1RMS = 0; VTerm1Phase = 0;
+		  ITerm2RMS = 1; ITerm2Phase = 0;
+	EndIf
+EndIf
+
+
+modelDim = 3; // 
+Flag_Axi = 0; // 1 for AXI - not tested for the moment, it makes sense only for modelDim = 2
+
+h2Ddepth = (modelDim==3) ? 1 : (Flag_Axi ? 2*Pi : l);
\ No newline at end of file
diff --git a/cylinderCoax/geo_pro/cylinder2Daxi.geo b/cylinderCoax/geo_pro/cylinder2Daxi.geo
new file mode 100644
index 0000000..1d3732e
--- /dev/null
+++ b/cylinderCoax/geo_pro/cylinder2Daxi.geo
@@ -0,0 +1,145 @@
+Include "cylinder2Daxi_data.pro";
+
+// hereafter there are two variants of code, both are correct
+
+delta = Sqrt(2.0/(2*Pi*Freq*mu*sigma));
+If (delta < a)
+	lcar_p0 = a - delta;
+    lcar_p1 = delta;
+	If (lcar_p0 < lcar_p1)
+		lcar_p1 = lcar_p0;
+	    Printf(" ------------------- > lcar_p0 < l_car_p1 at Freg = %g",Freq);
+	Else
+		Printf(" ==== > lcar_p0 > l_car_p1 at Freg = %g",Freq);
+	EndIf
+Else
+	Printf(" :::::::::::: delta > a at Freg = %g",Freq);
+	lcar_p0 = a;
+    lcar_p1 = a;
+EndIf
+
+lcar_p0 = s*lcar_p0 ;
+lcar_p1 = s*lcar_p1 ;
+
+If (_use_uniformMesh==1)
+    lcar = s*a/10;
+
+	pp[]+=newp; Point(newp) = {0, 0, 0, lcar};
+	pp[]+=newp; Point(newp) = {a, 0, 0, lcar};
+	pp[]+=newp; Point(newp) = {a, l, 0, lcar};
+	pp[]+=newp; Point(newp) = {0, l, 0, lcar};
+	For k In {0:3}
+		If (k == 3)
+			ll[]+=newl; Line(newl) = {pp[k],pp[0]};
+		Else
+			ll[]+=newl; Line(newl) = {pp[k],pp[k+1]};
+		EndIf
+	EndFor
+
+	cl = newll;  Line Loop(cl) = ll[];
+	surfRectangle = news; Plane Surface(surfRectangle) = {cl};
+	
+	If(_use_transfinite)
+		nopa = Ceil[a/lcar] ;
+		Transfinite Line {ll[0],ll[2]} = nopa; 
+		nopl = Ceil[l/lcar] ;
+		Transfinite Line {ll[1],ll[3]} = nopl; 
+		Transfinite Surface(surfRectangle);
+	EndIf
+
+	If (_use_recombine)
+		Recombine Surface{surfRectangle};
+	EndIf
+	
+	/* Definition of Physical entities (surfaces, lines). The Physical
+       entities tell GMSH the elements and their associated region numbers
+       to save in the file 'Ishape2d.msh'. */
+
+	Physical Surface ("MaterialX", MATERIALX) = {surfRectangle};     
+
+	Physical Line ("Ground",   GROUND) = {ll[2]};                 
+	Physical Line ("Terminal", TERMINAL) = {ll[0]};               
+	Physical Line ("RightBoundary", RIGHTBOUNDARY) = {ll[1]} ;     
+	Physical Line ("LeftBoundary", LEFTBOUNDARY) = {ll[3]} ;      
+
+	Physical Line ("Boundary", BOUNDARY) = {ll[]} ; 
+	
+Else
+
+    // non-uniform mesh
+	
+	pp[]+=newp; Point(newp) = {0, 0, 0, lcar_p0};
+	pp[]+=newp; Point(newp) = {a, 0, 0, lcar_p1};
+	pp[]+=newp; Point(newp) = {a, l, 0, lcar_p1};
+	pp[]+=newp; Point(newp) = {0, l, 0, lcar_p0};
+	For k In {0:3}
+		If (k == 3)
+			ll[]+=newl; Line(newl) = {pp[k],pp[0]};
+		Else
+			ll[]+=newl; Line(newl) = {pp[k],pp[k+1]};
+		EndIf
+	EndFor
+
+	cl = newll;  Line Loop(cl) = ll[];
+	surfRectangle = news; Plane Surface(surfRectangle) = {cl};
+	
+	Printf("================================ _use_progression = %g",_use_progression);
+	
+	
+    If (_use_progression == 1)
+		Printf("================================ delta = %g",delta);
+		    
+		pr1 =  0.5;
+		pr2 =  2;
+		
+		noptest = Log[1-a*(1-pr2)/delta]/Log[pr2]-1;
+		noptest = Ceil[noptest];
+		Printf("================================ noptest = %g",noptest);
+	    If (noptest < 4)
+					noptest = 4;
+		EndIf
+		Printf("================================ noptest2 = %g",noptest);
+ 		nopp = noptest;
+		Transfinite Line {ll[0]}= nopp*s Using Progression pr1; 
+		Transfinite Line {ll[2]}= nopp*s Using Progression pr2; 
+	EndIf	
+		
+		
+	
+	If(_use_transfinite)
+		nopl = 2; 
+		Transfinite Line {ll[1],ll[3]} = nopl; 
+		Transfinite Surface(surfRectangle);
+		If (_use_recombine)
+			Recombine Surface{surfRectangle};
+		EndIf
+	EndIf
+
+	/* Definition of Physical entities (surfaces, lines). The Physical
+       entities tell GMSH the elements and their associated region numbers
+       to save in the file 'Ishape2d.msh'. */
+
+	Physical Surface ("MaterialX", MATERIALX) = {surfRectangle};      
+
+	Physical Line ("Ground",   GROUND) = {ll[2]};                 
+	Physical Line ("Terminal", TERMINAL) = {ll[0]};               
+	Physical Line ("RightBoundary", RIGHTBOUNDARY) = {ll[1]} ;     
+	Physical Line ("LeftBoundary", LEFTBOUNDARY) = {ll[3]} ;       
+
+	Physical Line ("Boundary", BOUNDARY) = {ll[]} ; 
+	
+EndIf
+
+If(_use_transfinite)
+	If (_use_recombine)
+		Recursive Color Red { Surface{surfRectangle};}
+	Else
+		Recursive Color Blue { Surface{surfRectangle};}
+	EndIf
+Else
+	If (_use_recombine)
+		Recursive Color Orange { Surface{surfRectangle};}
+	Else
+		Recursive Color Black { Surface{surfRectangle};}
+	EndIf
+EndIf
diff --git a/cylinderCoax/geo_pro/cylinder2Daxi.pro b/cylinderCoax/geo_pro/cylinder2Daxi.pro
new file mode 100644
index 0000000..4832b2a
--- /dev/null
+++ b/cylinderCoax/geo_pro/cylinder2Daxi.pro
@@ -0,0 +1,175 @@
+Include "cylinder2Daxi_data.pro";
+
+Group{
+   MaterialX = Region[MATERIALX]; /* MaterialX - homogenous domain */
+   Ground = Region[GROUND];       /* Boundary - various parts */
+   Terminal = Region[TERMINAL];
+   RightBoundary = Region[RIGHTBOUNDARY];
+   Axis = Region[LEFTBOUNDARY];                 // the axis is taken out of the domain
+   Cut1H = Region[RIGHTBOUNDARY];               // for the H formulation
+   BoundaryNotTerminal = Region[{RIGHTBOUNDARY}];
+  
+   Sur_Terminals_FWece = Region[{Ground, Terminal}]; // all terminals
+  
+   Vol_FW = Region[{MaterialX}]; // domain without the boundary
+   Sur_FW = Region[{Sur_Terminals_FWece, RightBoundary}]; //all boundary
+   Dom_FW = Region[ {Vol_FW, Sur_FW, Axis} ];
+ }
+ 
+ Function{
+  nu[#{MaterialX,Sur_FW}] = nu;             // magnetic reluctivity [m/H]
+  mu[#{MaterialX,Sur_FW}] = mu;             // permeability [H/m]
+  epsilon[#{MaterialX,Sur_FW}] = eps;       // permittivity [F/m]
+  sigma[#{MaterialX,Sur_FW}]  = sigma;      // conductivity [S/m]
+
+  // This is a problem with 2 terminals, out of which one is grounded.
+  // The non-grounded can be excited in voltage or in current.
+
+  If (Flag_AnalysisType == 0)    // voltage excitation 
+    // (Flag_AnalysisTypeFormulation==0) => formulation in E
+    VTerminal1[] = (Flag_AnalysisTypeFormulation==0) ? -1 : 1;
+  Else                            // current excitation
+    ITerminal1[] = (Flag_AnalysisTypeFormulation==0) ? -1 : 1;
+  EndIf
+
+ }
+
+Constraint{
+
+  // ece BC
+  { Name SetTerminalPotential; Type Assign;  // voltage excited terminals
+    Case {
+      { Region Ground; Value 0.; }
+      If((Flag_AnalysisType==0))    
+        { Region Terminal; Value VTerminal1[]; }
+      EndIf
+    }
+  }
+  { Name SetTerminalCurrent; Type Assign;   // current excited terminals
+    Case {
+      If((Flag_AnalysisType==1))
+        { Region Terminal; Value ITerminal1[]/h2Ddepth; }  // here the depth is needed
+      EndIf
+    }
+  }
+  { Name SetTerminalCurrentH; Type Assign;   // current excited terminals - formulation in H
+    Case {
+      If(Flag_AnalysisType==1)
+        { Region Cut1H; Value ITerminal1[]/h2Ddepth; }  
+      EndIf
+    }
+  }
+   
+  { Name SetTerminalPotentialH ;
+    Case {
+    If(Flag_AnalysisType==0)
+      { Region Cut1H; Value VTerminal1[] ; }
+    EndIf
+    }
+  }
+  
+  { Name ImposeHonAxis ;
+    Case {
+      If( Flag_AnalysisTypeFormulation==1 && Flag_Axi==1 ) // H
+          { Region Axis; Value 0 ; }
+      EndIf
+    }
+	}
+	
+  { Name ImposeHorder2onAxis ;
+    Case {
+      If( Flag_AnalysisTypeFormulation==1 && Flag_Axi==1 ) // H
+        { Region Axis; Value 0 ; }
+      EndIf
+    }
+  }
+   
+  { Name ImposeE ;
+  }
+
+}
+
+modelDir = "../";
+Dir = "res/";
+
+// The formulation and its tools
+If (Flag_AnalysisTypeFormulation==0)          
+  Include "../formulations/FullWave_E_ece_SISO_secondOrder.pro"                    // formulation in E
+Else 
+  Include "../formulations/FullWave_H_ece_SISO_cplx_2D_and_AXI_secondOrder.pro"    // formulation in H
+EndIf
+
+ 
+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
+
+PostOperation {
+  If (Flag_AnalysisTypeFormulation==0)   // E
+
+    { Name Maps; NameOfPostProcessing FW_E_ece ;
+      Operation {
+        Print [ gradV, OnElementsOf ElementsOf[Vol_FW,OnOneSideOf Sur_FW], File StrCat[Dir,"map_gradV.pos" ]];
+        Print [ Vterminal, 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 [ B, OnElementsOf Vol_FW, File StrCat[Dir,"map_B.pos" ]];
+        Print [ H, OnElementsOf Vol_FW, File StrCat[Dir,"map_H.pos" ]];
+
+        Print [ VnotTerminals, OnElementsOf Vol_FW, File StrCat[Dir,"map_VnotTerminals.pos" ]];
+
+        Call Change_post_options;
+      }
+    }   
+
+    { Name TransferMatrix; NameOfPostProcessing FW_E_ece ;
+      Operation {
+        If(Flag_AnalysisType==1)  // current excitation
+          Print [ Vterminal, OnRegion Terminal,
+            Format Table , File  > StrCat[Dir, "test_Z_RI_formE.s1p"]] ;
+        Else
+          Print [ I, OnRegion Terminal,
+            Format Table , File  > StrCat[Dir, "test_Y_RI_formE.s1p"]] ;
+        EndIf
+      }  
+    }
+  
+  Else // H
+
+    { Name Maps; NameOfPostProcessing FW_H_ece ;
+      Operation {
+        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 [ B, OnElementsOf Vol_FW, File StrCat[Dir,"map_B.pos" ]];
+        Print [ H, OnElementsOf Vol_FW, File StrCat[Dir,"map_H.pos" ]];
+        Call Change_post_options;
+      }    
+    }
+
+    { Name TransferMatrix; NameOfPostProcessing FW_H_ece ;
+      Operation {
+        If(Flag_AnalysisType==1)  // current excitation
+          Print [ Vterminal, OnRegion Cut1H,
+            Format Table , File  > StrCat[Dir, "test_Z_RI_formH.s1p"]] ;
+        Else
+          Print [ I, OnRegion Cut1H,
+             Format Table , File  > StrCat[Dir, "test_Y_RI_formH.s1p"]] ;
+        EndIf
+      }  
+    }  
+
+  EndIf
+}
+
diff --git a/cylinderCoax/geo_pro/cylinder2Daxi_data.pro b/cylinderCoax/geo_pro/cylinder2Daxi_data.pro
new file mode 100644
index 0000000..0e10b24
--- /dev/null
+++ b/cylinderCoax/geo_pro/cylinder2Daxi_data.pro
@@ -0,0 +1,202 @@
+/*
+Geometry:
+a = radius of the cylinder [m]
+l = height of the cylinder [m]
+==
+Material 
+epsr = relative permittivity
+mur = relative permeability
+sigma = conductivity [S/m]
+*/
+
+cGUI= "{TEST-Cilinder2Daxi-ProblemData/";
+cGUI0 = StrCat[cGUI, "0"];
+cGUI1 = StrCat[cGUI, "1"];
+cGUI2 = StrCat[cGUI, "2"];
+cGUI3 = StrCat[cGUI, "3"];
+cGUI4 = StrCat[cGUI, "4"];
+close_menu = 1;
+
+/* new menu for the interface - which test */
+mTypeTest = StrCat[cGUI0,"Test/0"];
+colorMTypeTest = "Orange";
+/* new menu for the interface - frequencies */
+mValuesBC = StrCat[cGUI1,"Frequency/0"];
+colorMValuesBC = "Ivory";
+/* new menu for the interface - type of terminal excitation */
+mTypeBC = StrCat[cGUI2,"Type of BC/0"];
+colorMTypeBC = "Red";
+/* new menu for the interface - type of formulation */
+mTypeBC2 = StrCat[cGUI3,"Type of formulation/0"];
+colorMTypeBC = "Blue";
+/* new menu for the interface - type of formulation */
+mTypeMesh = StrCat[cGUI4,"MeshSettings/"];
+colorMTypeMesh = "Red";
+
+/* Test */
+DefineConstant[
+  Flag_WhichTest = { 0,
+    Choices{
+      0="a = 2.5 um, l = 10 um",
+      1="a = 4 mm, l = 10 mm"},
+    Name Str[mTypeTest], Highlight Str[colorMTypeTest], Closed !close_menu,
+	ServerAction Str["Reset", "GetDP/1ResolutionChoices" , ',' , StrCat[mValuesBC, "0Working Frequency"] ]}
+	];
+	
+If (Flag_WhichTest==0)
+	/* dimensions used in the SCEE 2020 and JMI papers */
+	a = 2.5e-6;
+	l = 10e-6;
+	fmin = 1e7;   // Hz
+	fmax = 100e9; // Hz
+	nop = 20;
+    freqs()=LogSpace[Log10[fmin],Log10[fmax],nop];
+    //freqs()=LinSpace[fmin,fmax,nop];
+Else
+	/* a inspired  from Wu coax, but l is short */
+	a = 4e-3;
+	l = 10e-3;
+	fmin = 1;   // Hz
+	//fmax = 600e6; // Hz
+	fmax = 100e6; // Hz
+	//fmax = 1e4; 
+	nop = 20;
+	freqs()=LogSpace[Log10[fmin],Log10[fmax],nop];
+    //freqs()=LinSpace[fmin,fmax,nop];
+EndIf
+	
+epsr = 1;
+mur = 1;
+sigma = 6.6e7;
+
+eps0 = 8.854187818e-12; // F/m - absolute permitivity of vacuum
+mu0  = 4.e-7 * Pi;      // H/m - absolute permeability of vacuum
+nu0  = 1/mu0;           // m/H - absolute reluctivity of vacuum
+
+eps = eps0*epsr;        // F/m - absolute permitivity
+mu = mu0*mur;           // F/m - absolute permeability
+nu = 1/mu;   	
+	
+/* Frequencies */
+DefineConstant[
+_use_freq_loop = {0,  Choices{0,1}, Name StrCat[mValuesBC, "0Loop on frequencies?"],
+    ServerAction Str["Reset", StrCat[mValuesBC, "0Working Frequency"]], Highlight "Yellow"}
+  Freq  = {freqs(0), Choices{freqs()}, Loop _use_freq_loop, Name StrCat[mValuesBC, "0Working Frequency"],
+    Units "Hz", Highlight  Str[colorMValuesBC], Closed !close_menu}
+];
+
+/* Excitation type */
+DefineConstant[
+  Flag_AnalysisType = {  TERMINAL_EXC,
+    Choices{
+      0="bottomTerm ev",
+      1="bottomTerm ec"},
+    Name Str[mTypeBC], Highlight Str[colorMTypeBC], Closed !close_menu,
+    ServerAction Str["Reset", "GetDP/1ResolutionChoices"]}
+	];
+	
+	
+	
+
+/* Which formulation */
+DefineConstant[
+  Flag_AnalysisTypeFormulation = { FORMULATION,
+    Choices{
+      0="0: formulation in E inside + ECE",
+      1="1: formulation in H inside + ECE"
+    },
+    Name Str[mTypeBC2], Highlight Str[colorMTypeBC], Closed !close_menu,
+    Help Str["E inside, electric V on the boundary",
+      "H inside, reduced magnetic V on the boundary"],
+    ServerAction Str["Reset", "GetDP/1ResolutionChoices"]}
+];
+
+//s_min = 0.5; // the finest mesh
+s_min = 1;
+//s_min = 0.7071067811865475; // 1/sqrt(2)
+//s_min = 1.414213562373095; //sqrt(2)
+//
+//s_min = 7.75E-3;
+//s_min = 3.31E-2;
+
+s_max = 1.414213562373095; //the coarsest mesh
+nops = 7;
+
+
+s_values()=LinSpace[s_min,s_max,nops];
+
+
+/* Mesh setings */
+DefineConstant[
+   _use_uniformMesh = {1, Choices{0,1}, Name StrCat[mTypeMesh,"0Uniform mesh?"], Highlight "Pink", Closed !close_menu},
+   _use_progression = {0, Choices{0,1}, Name StrCat[mTypeMesh,"1Use geometric progression on the bottom line (only for non-uniform mesh)?"], Highlight "Pink", Closed !close_menu, ReadOnly _use_uniformMesh},
+   _use_transfinite = {0, Choices{0,1}, Name StrCat[mTypeMesh,"2One element along Oy ; Transfinite for uniform mesh?"], Highlight "Pink", Closed !close_menu, ReadOnly _use_uniformMesh},
+   _use_recombine = {0, Choices{0,1}, Name StrCat[mTypeMesh,"3Use recombine"], Highlight "Pink", Closed !close_menu},
+   //_use_s_loop = {0, Choices{0,1}, Name StrCat[mTypeMesh,"4Use s loop"],
+   // ServerAction Str["Reset", StrCat[mTypeMesh, "5Mesh factor s, where car.length = s*a*0.1"]],   Highlight "Pink", Closed !close_menu},
+   s = {s_min, Name StrCat[mTypeMesh,"6Mesh factor"]}
+   //s  = {s_values(0), Choices{s_values()}, Loop _use_s_loop, Name StrCat[mTypeMesh, "6Mesh factor s, where car.length = s*a*0.1"],
+   // Units "-", Highlight  Str[colorMTypeMesh], Closed !close_menu}
+   FE_Order = { FE_ORDER,
+    Choices{
+      1="1: first order elements",
+      2="2: second order elements"
+    },
+    Name StrCat[mTypeMesh,"7Element order"], Highlight Str[colorMTypeMesh], Closed !close_menu,
+    ServerAction Str["Reset", "GetDP/1ResolutionChoices"]}
+   
+];
+
+Printf("a = %g  m",a);
+Printf("l = %g  m",l);
+Printf("fmin = %g  Hz",fmin);
+Printf("fmax = %g  Hz",fmax);
+Printf("s = %g ",s);
+Printf("FE order = %g ",FE_Order);
+
+
+modelDim = 2; // 
+Flag_Axi = 1; // 1 for AXI - it makes sense only for modelDim = 2
+
+If ((modelDim == 2)&&(Flag_Axi == 0)) // 2D
+	h2Ddepth = h;
+ElseIf ((modelDim == 2)&&(Flag_Axi == 1))   // 2D AXI // 2.5 D
+	h2Ddepth = 2*Pi;
+Else // 3D
+    h2Ddepth = 1;
+EndIf
+
+//=============================
+// Indexes of physical entities
+// Surfaces
+MATERIALX = 100;
+
+// boundaries
+GROUND   = 120 ;
+TERMINAL = 121 ;
+RIGHTBOUNDARY = 131 ;
+LEFTBOUNDARY  = 132; // axis
+
+BOUNDARY  = 331; 
+
+
+
+/* Hint:
+
+Concerning the 2nd order elements, there are two ways of doing it:
+
+1) using hierarchical elements (BFs associated to nodes and edges; or to
+edges and facets);
+
+2) using 2nd order mesh and the corresponding BFs => adding in the geo
+file: Mesh.ElementOrder = 2; // Element order (1: first order elements)
+
+For the former, you need to modify add some lines in the pro-file of the BFs
+For the latter, you have to generate the 2nd order mesh and add the new
+second order elements in the integration.
+
+1) and 2) are incompatible => 1) works with 1st order geometrical elements
+2) it is experimental, element-type: Triangle2, Quadrangle2,
+Tetrahedron2 and so forth.
+
+*/
\ No newline at end of file
diff --git a/cylinderCoax/geo_pro/cylinder3D.geo b/cylinderCoax/geo_pro/cylinder3D.geo
new file mode 100644
index 0000000..78861ee
--- /dev/null
+++ b/cylinderCoax/geo_pro/cylinder3D.geo
@@ -0,0 +1,89 @@
+Include "cylinder3D_data.pro";
+
+delta = Sqrt(2.0/(2*Pi*Freq*mu*sigma));
+Printf("Freq = %e",Freq);
+If (delta < a)
+	Printf("delta less than a");
+	lcar = delta/nbDelta;
+Else
+	Printf("delta GREATER than a");
+	lcar = l/10;
+EndIf
+
+
+lcarCenter = a/3;
+
+p0 = newp; Point(p0) = {0,0,0,lcarCenter};
+p1 = newp; Point(p1) = {a,0,0,lcar};
+p2 = newp; Point(p2) = {0,a,0,lcar};
+p3 = newp; Point(p3) = {-a,0,0,lcar};
+p4 = newp; Point(p4) = {0,-a,0,lcar};
+
+p0l = newp; Point(p0l) = {0,0,l,lcarCenter};
+p1l = newp; Point(p1l) = {a,0,l,lcar};
+p2l = newp; Point(p2l) = {0,a,l,lcar};
+p3l = newp; Point(p3l) = {-a,0,l,lcar};
+p4l = newp; Point(p4l) = {0,-a,l,lcar};
+
+
+c1 = newc; Circle(c1) = {p1,p0,p2};
+c2 = newc; Circle(c2) = {p2,p0,p3};
+c3 = newc; Circle(c3) = {p3,p0,p4};
+c4 = newc; Circle(c4) = {p4,p0,p1};
+
+c1l = newc; Circle(c1l) = {p1l,p0l,p2l};
+c2l = newc; Circle(c2l) = {p2l,p0l,p3l};
+c3l = newc; Circle(c3l) = {p3l,p0l,p4l};
+c4l = newc; Circle(c4l) = {p4l,p0l,p1l};
+
+l1 = newl; Line(l1) = {p1, p1l};
+l2 = newl; Line(l2) = {p2, p2l};
+l3 = newl; Line(l3) = {p3, p3l};
+l4 = newl; Line(l4) = {p4, p4l};
+
+cL = newll;   Curve Loop(cL) = {c1,c2,c3,c4};      scL = news;  Surface(scL) = {cL};
+cLl = newll;  Curve Loop(cLl) = {c1l,c2l,c3l,c4l}; scLl = news; Surface(scLl) = {cLl};
+cL12 = newll; Curve Loop(cL12) = {c1,l2,-c1l,-l1}; scL12 = news; Surface(scL12) = {cL12};
+cL23 = newll; Curve Loop(cL23) = {c2,l3,-c2l,-l2}; scL23 = news; Surface(scL23) = {cL23};
+cL34 = newll; Curve Loop(cL34) = {c3,l4,-c3l,-l3}; scL34 = news; Surface(scL34) = {cL34};
+cL41 = newll; Curve Loop(cL41) = {c4,l1,-c4l,-l4}; scL41 = news; Surface(scL41) = {cL41};
+
+closedS = newsl; Surface Loop(closedS) = {scL,scLl,scL12,scL23,scL34,scL41};
+volDom = newv; Volume(volDom) = {closedS};
+
+Transfinite Line{l1} = 10*s;
+Transfinite Line{l2} = 10*s;
+Transfinite Line{l3} = 10*s;
+Transfinite Line{l4} = 10*s;
+Transfinite Surface {scL12};
+Transfinite Surface {scL23};
+Transfinite Surface {scL34};
+Transfinite Surface {scL41};
+
+Field[1] = Distance;
+Field[1].SurfacesList = {scL12,scL23,scL34,scL41};
+Field[1].Sampling = 50;
+Field[2] = Threshold;
+Field[2].InField = 1;
+Field[2].SizeMin = lcar;
+Field[2].SizeMax = a/3;
+Field[2].DistMin = 0;
+Field[2].DistMax = 3*delta;
+
+Background Field = 2; 
+Mesh.MeshSizeExtendFromBoundary = 0;
+Mesh.MeshSizeFromPoints = 0;
+Mesh.MeshSizeFromCurvature = 0;
+Mesh.Algorithm = 5;
+
+
+
+Physical Volume ("MaterialX", 100) = {volDom} ;      
+
+Physical Surface ("Ground",   120) = {scLl} ;          
+Physical Surface ("Terminal", 121) = {scL} ;          
+Physical Surface ("BoundaryNotTerminal", 131) = {scL12,scL23,scL34,scL41} ;  
+
+// cut for the H formulation
+Cohomology(1) {{131},{}}; // bnd fara term
+CUT1H = 132;
\ No newline at end of file
diff --git a/cylinderCoax/geo_pro/cylinder3D.pro b/cylinderCoax/geo_pro/cylinder3D.pro
new file mode 100644
index 0000000..1822241
--- /dev/null
+++ b/cylinderCoax/geo_pro/cylinder3D.pro
@@ -0,0 +1,160 @@
+Include "cylinder3D_data.pro";
+
+Group{
+   MaterialX = Region[MATERIALX]; /* MaterialX - homogenous domain */
+   Ground = Region[GROUND];       /* Boundary - various parts */
+   Terminal = Region[TERMINAL];
+   
+   BoundaryNotTerminal = Region[BOUNDARYNOTTERMINAL];
+   Cut1H = Region[132];  // for the H formulation  // this is CUT1H
+
+   Vol_FW = Region[{MaterialX}]; // domain without the boundary
+   Sur_Terminals_FWece = Region[{Ground, Terminal}]; // all terminals
+   Sur_FW = Region[{Sur_Terminals_FWece, BoundaryNotTerminal}]; //all boundary
+   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 Ishape2d_data.pro */
+
+   nu[#{MaterialX,Sur_FW}] = nu;  // magnetic reluctivity [m/H]
+   mu[#{MaterialX,Sur_FW}] = mu;  // permeability [H/m]
+   epsilon[#{MaterialX,Sur_FW}] = eps;       // permittivity [F/m]
+   sigma[#{MaterialX,Sur_FW}]  = sigma;      // conductivity [S/m]
+
+   // FW in E, ece BC, voltage excitation on bottom, ground on top
+   // When imposing voltage or current via a constraint, the region appears there
+   // If you have a complex, the value  has to be a function => use square brackets []
+   If (Flag_AnalysisTypeFormulation==0) // E
+      If(Flag_AnalysisType==0) // voltage excitation
+         VTerminal1[] = -VTermRMS*Complex[Cos[VTermPhase],Sin[VTermPhase]];
+      Else  // current excitation
+         ITerminal1[] = -ITermRMS*Complex[Cos[ITermPhase],Sin[ITermPhase]];
+      EndIf
+   Else // H
+      If(Flag_AnalysisType==0) // voltage excitation
+         VTerminal1[] = VTermRMS*Complex[Cos[VTermPhase],Sin[VTermPhase]];
+      Else  // current excitation
+         ITerminal1[] = ITermRMS*Complex[Cos[ITermPhase],Sin[ITermPhase]];
+      EndIf
+   EndIf
+ }
+
+ Constraint{
+
+   // ece BC or electrokinetics
+   { Name SetTerminalPotential; Type Assign;  // voltage excited terminals
+     Case {
+       { Region Ground; Value 0.; }
+       If(Flag_AnalysisType==0)
+         { Region Terminal; Value VTerminal1[]; }
+       EndIf
+     }
+   }
+   { Name SetTerminalCurrent; Type Assign;   // current excited terminals
+     Case {
+       If(Flag_AnalysisType==1)
+         { Region Terminal; Value ITerminal1[]/h2Ddepth; }  
+       EndIf
+     }
+   }
+   
+   { Name SetTerminalCurrentH; Type Assign;   // current excited terminals - formulation in H
+     Case {
+       If(Flag_AnalysisType==1)
+         { Region Cut1H; Value ITerminal1[]/h2Ddepth; }  
+       EndIf
+     }
+   }
+   
+   { Name SetTerminalPotentialH ;
+    Case {
+      If(Flag_AnalysisType==0)
+      { Region Cut1H; Value VTerminal1[] ; }
+      EndIf
+    }
+  }
+
+}
+
+modelDir = "../";
+Dir = "res/";
+
+// The formulation and its tools
+If (Flag_AnalysisTypeFormulation==0)  // E inside
+  Include "../formulations/FullWave_E_ece_SISO_secondOrder.pro"
+Else // H inside
+  Include "../formulations/FullWave_H_ece_SISO_cplx_3D_secondOrder.pro"
+EndIf
+
+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
+
+
+PostOperation {
+  If (Flag_AnalysisTypeFormulation==0)   // E
+
+    { Name Maps; NameOfPostProcessing FW_E_ece ;
+      Operation {
+        Print [ gradV, OnElementsOf ElementsOf[Vol_FW,OnOneSideOf Sur_FW], File StrCat[Dir,"map_gradV.pos" ]];
+        Print [ Vterminal, 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 [ B, OnElementsOf Vol_FW, File StrCat[Dir,"map_B.pos" ]];
+        Print [ H, OnElementsOf Vol_FW, File StrCat[Dir,"map_H.pos" ]];
+        Print [ VnotTerminals, OnElementsOf Vol_FW, File StrCat[Dir,"map_VnotTerminals.pos" ]];
+
+        Call Change_post_options;
+      }
+    }
+
+    { Name TransferMatrix; NameOfPostProcessing FW_E_ece ;
+      Operation {
+        If(Flag_AnalysisType==1)  // current excitation        
+          Print [ Vterminal, OnRegion Terminal,
+            Format Table , File  > StrCat[Dir, "test_Z_RI_formE.s1p"]] ;
+        Else 
+          Print [ I, OnRegion Terminal,
+            Format Table , File  > StrCat[Dir, "test_Y_RI_formE.s1p"]] ;
+        EndIf
+      }
+    }
+  
+  Else // H
+ 
+    { Name TransferMatrix; NameOfPostProcessing FW_H_ece ;
+      Operation {
+        If(Flag_AnalysisType==1)  // current excitation        
+          Print [ Vterminal, OnRegion Cut1H,
+            Format Table , File  > StrCat[Dir, "test_Z_RI_formH.s1p"]] ;
+          Else      
+          Print [ I, OnRegion Cut1H,
+            Format Table , File  > StrCat[Dir, "test_Y_RI_formH.s1p"]] ;
+        EndIf
+      }
+    }
+
+    { Name Maps; NameOfPostProcessing FW_H_ece ;
+      Operation {
+        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 [ B, OnElementsOf Vol_FW, File StrCat[Dir,"map_B.pos" ]];
+        Print [ H, OnElementsOf Vol_FW, File StrCat[Dir,"map_H.pos" ]];
+        Call Change_post_options;
+      }
+    }  
+  EndIf
+
+}
\ No newline at end of file
diff --git a/cylinderCoax/geo_pro/cylinder3D_data.pro b/cylinderCoax/geo_pro/cylinder3D_data.pro
new file mode 100644
index 0000000..310098b
--- /dev/null
+++ b/cylinderCoax/geo_pro/cylinder3D_data.pro
@@ -0,0 +1,167 @@
+/*
+Geometry:
+a = radius of the cylinder [m]
+l = height of the cylinder [m]
+==
+Material 
+epsr = relative permittivity
+mur = relative permeability
+sigma = conductivity [S/m]
+*/
+colorro = "LightGrey";
+
+cGUI= "{TEST-Cilinder3D-ProblemData/";
+cGUI0 = StrCat[cGUI, "0"];
+cGUI1 = StrCat[cGUI, "1"];
+cGUI2 = StrCat[cGUI, "2"];
+cGUI3 = StrCat[cGUI, "3"];
+cGUI4 = StrCat[cGUI, "4"];
+close_menu = 1;
+
+/* new menu for the interface - which test */
+mTypeTest = StrCat[cGUI0,"Test/0"];
+colorMTypeTest = "Orange";
+/* new menu for the interface - frequencies */
+mValuesBC = StrCat[cGUI1,"Frequency/0"];
+colorMValuesBC = "Ivory";
+/* new menu for the interface - type of terminal excitation */
+mTypeBC = StrCat[cGUI2,"Type of BC/0"];
+colorMTypeBC = "Red";
+/* new menu for the interface - type of formulation */
+mTypeBC2 = StrCat[cGUI3,"Type of formulation/0"];
+colorMTypeBC = "Blue";
+/* new menu for the interface - type of formulation */
+mTypeMesh = StrCat[cGUI4,"MeshSettings/"];
+colorMTypeMesh = "Red";
+
+/* Test */
+DefineConstant[
+  Flag_WhichTest = { 0,
+    Choices{
+      0="a = 2.5 um, l = 10 um",
+      1="a = 4 mm, l = 10 mm"},
+    Name Str[mTypeTest], Highlight Str[colorMTypeTest], Closed !close_menu,
+	ServerAction Str["Reset", "GetDP/1ResolutionChoices" , ',' , StrCat[mValuesBC, "0Working Frequency"] ]}
+	];
+	
+If (Flag_WhichTest==0)
+	/* dimensions used in the SCEE 2020, 2022 and JMI papers */
+	a = 2.5e-6;
+	l = 10e-6;
+	fmin = 1e7;   // Hz
+	fmax = 100e9; // Hz
+	nop = 10;
+    //freqs()=LogSpace[Log10[fmin],Log10[fmax],nop];
+    freqs()=LinSpace[fmin,fmax,nop];
+Else
+	/* a inspired  from Wu coax, but l is short */
+	a = 4e-3;
+	l = 10e-3;
+	fmin = 1;   // Hz
+	//fmax = 600e6; // Hz
+	fmax = 100e6; // Hz
+	//fmax = 1e4; 
+	nop = 20;
+	freqs()=LogSpace[Log10[fmin],Log10[fmax],nop];
+    //freqs()=LinSpace[fmin,fmax,nop];
+EndIf
+
+epsr = 1;
+mur = 1;
+sigma = 6.6e7;
+
+eps0 = 8.854187818e-12; // F/m - absolute permitivity of vacuum
+mu0  = 4.e-7 * Pi;      // H/m - absolute permeability of vacuum
+nu0  = 1/mu0;           // m/H - absolute reluctivity of vacuum
+
+eps = eps0*epsr;        // F/m - absolute permitivity
+mu = mu0*mur;           // F/m - absolute permeability
+nu = 1/mu;   	
+	
+// Frequencies
+DefineConstant[
+_use_freq_loop = {0,  Choices{0,1}, Name StrCat[mValuesBC, "0Loop on frequencies?"],
+    ServerAction Str["Reset", StrCat[mValuesBC, "0Working Frequency"]], Highlight "Yellow"}
+  Freq  = {freqs(0), Choices{freqs()}, Loop _use_freq_loop, Name StrCat[mValuesBC, "0Working Frequency"],
+    Units "Hz", Highlight  Str[colorMValuesBC], Closed !close_menu}
+];
+
+// Excitation type
+DefineConstant[
+  Flag_AnalysisType = { TERMINAL_EXC,
+    Choices{
+      0="bottomTerm ev",
+      1="bottomTerm ec"},
+    Name Str[mTypeBC], Highlight Str[colorMTypeBC], Closed !close_menu,
+    ServerAction Str["Reset", "GetDP/1ResolutionChoices"]}
+	];
+
+
+// Which formulation
+DefineConstant[
+  Flag_AnalysisTypeFormulation = { FORMULATION,
+    Choices{
+      0="0: formulation in E inside + ECE",
+      1="1: formulation in H inside + ECE"
+    },
+    Name Str[mTypeBC2], Highlight Str[colorMTypeBC], Closed !close_menu,
+    Help Str["E inside, electric V on the boundary",
+      "H inside, reduced magnetic V on the boundary"],
+    ServerAction Str["Reset", "GetDP/1ResolutionChoices"]}
+];
+
+
+meshSettings = StrCat[cGUI,"5MeshSettings/0"];
+
+DefineConstant[
+   s = {1, Name StrCat[meshSettings, "1Mesh factor/1Mesh factor for Oz (1 means 10 nodes on Oz)"]}
+   nbDelta = {1, Name StrCat[meshSettings, "1Mesh factor/2Nb of elems per skin depth (lcar is delta over nbDelta)"]}
+   Flag_ElementOrder = { FE_ORDER,
+    Choices{
+      1="1: first order elements",
+      2="2: second order elements"
+    },
+    Name StrCat[meshSettings, "2Element order"], Highlight "Red", Closed !close_menu,
+    ServerAction Str["Reset", "GetDP/1ResolutionChoices"]}
+	_use_recombine = {0, Choices{0,1}, Name StrCat[meshSettings, "3Use recombine?"], Highlight "Yellow", Closed !close_menu}
+];
+
+FE_Order = Flag_ElementOrder;
+modelDim = 3; 
+Flag_Axi = 0; // 1 for AXI - it makes sense only for modelDim = 2
+
+h2Ddepth = (modelDim==3) ? 1 : (Flag_Axi ? 2*Pi : l);
+
+DefineConstant[
+  // Flag_AnalysisType
+  // ==0 bottom ev
+  // ==1 bottom ec
+
+  VGroundRMS  = {0, Name StrCat[mValuesBC, "2Top/1Potential on bottom (rms)"], Units "V" , ReadOnly 1, Highlight Str[colorro],
+    Visible (Flag_AnalysisType==0) || (Flag_AnalysisType==1)}
+  VGroundPhase  = {0, Name StrCat[mValuesBC, "2Top/1Potential on bottom, phase"], Units "rad",  ReadOnly 1, Highlight Str[colorro],
+    Visible (Flag_AnalysisType==0) || (Flag_AnalysisType==1)}
+	
+  VTermRMS  = {1, Name StrCat[mValuesBC, "1Bottom/1Potential on Term1 (rms)"], Units "V" , ReadOnly 0, Highlight Str[colorMValuesBC],
+    Visible (Flag_AnalysisType==0)}
+  VTermPhase  = {0, Name StrCat[mValuesBC, "1Bottom/1Potential on Term1, phase"], Units "rad",  ReadOnly 0, Highlight Str[colorMValuesBC],
+    Visible (Flag_AnalysisType==0)}
+	
+  ITermRMS  = {1, Name StrCat[mValuesBC, "1Bottom/1Current through Term1 (enters the domain) (rms)"], Units "A" , ReadOnly 0, Highlight Str[colorMValuesBC],
+    Visible (Flag_AnalysisType==1) }
+  ITermPhase  = {0, Name StrCat[mValuesBC, "1Bottom/Current through Term1 (enters the domain), phase"], Units "rad",  ReadOnly 0, Highlight Str[colorMValuesBC],
+    Visible (Flag_AnalysisType==1) } 		
+];
+
+
+//=============================
+// Indexes of physical entities
+// Surfaces
+MATERIALX = 100;
+
+// boundaries
+GROUND   = 120 ;
+TERMINAL = 121 ;
+BOUNDARYNOTTERMINAL = 131 ;
+
+BOUNDARY  = 331; 
-- 
GitLab