diff --git a/ElectromagneticScattering/scattererTmatrix_data.geo b/ElectromagneticScattering/scattererTmatrix_data.geo
new file mode 100644
index 0000000000000000000000000000000000000000..8699e8983744a091acea15fe287c1a3a06666bfb
--- /dev/null
+++ b/ElectromagneticScattering/scattererTmatrix_data.geo
@@ -0,0 +1,177 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattererTmatrix_data.geo //
+///////////////////////////////
+
+nm = 1.;
+epsilon0 = 8.854187817e-3*nm;
+mu0 = 400.*Pi*nm;
+cel = 1.0/(Sqrt[epsilon0 * mu0]);
+deg2rad = Pi/180.;
+
+pp0 = "1Geometry/0";
+pp1 = "0Study Type/0";
+pp2 = "3Electromagnetic parameters/0";
+pp3 = "4Mesh size and PMLs parameters/0";
+pp4 = "5Postpro options/0";
+pp5 = "6Post plot options/0";
+close_menu = 0;
+colorro = "LightGrey";
+colorppOK = "Ivory";
+colorppWA = "LightSalmon1";
+colorppNO = "LightSalmon4";
+
+ELL = 0;
+PARALL = 1;
+CYL = 2;
+CONE = 3;
+TOR = 4;
+
+RES_PW = 0;
+RES_TMAT = 1;
+RES_GREEN = 2;
+RES_QNM = 3;
+
+DefineConstant[
+ flag_shape = {ELL , Name StrCat[pp0, "0Scatterer shape"],
+ Choices {ELL="ellispoid",PARALL="parallelepiped",CYL="cylinder",CONE="cone",TOR="split torus"},Closed 0}];
+
+If (flag_shape==ELL)
+ DefineConstant[
+ ell_rx = {73.45 , Name StrCat[pp0 , "1ellipsoid X-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ ell_ry = {73.45 , Name StrCat[pp0 , "2ellipsoid Y-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ ell_rz = {73.45 , Name StrCat[pp0 , "3ellipsoid Z-radius [nm]"], Highlight Str[colorppOK] , Closed 0}
+ ];
+ // FIXME gmsh Max?
+ If (ell_rx>ell_ry)
+ rbb_tmp=ell_rx;
+ Else rbb_tmp=ell_ry;
+ EndIf
+ If (ell_rz>rbb_tmp)
+ rbb=ell_rz;
+ Else rbb=rbb_tmp;
+ EndIf
+EndIf
+If (flag_shape==PARALL)
+ DefineConstant[
+ par_ax = {300 , Name StrCat[pp0 , "1cube X-edge size [nm]"], Highlight Str[colorppOK] , Closed 0},
+ par_ay = {100 , Name StrCat[pp0 , "2cube Y-edge size [nm]"], Highlight Str[colorppOK] , Closed 0},
+ par_az = {600 , Name StrCat[pp0 , "3cube Z-edge size [nm]"], Highlight Str[colorppOK] , Closed 0}];
+ rbb = 0.5*Sqrt[par_ax^2+par_ay^2+par_az^2];
+EndIf
+If (flag_shape==CYL)
+ DefineConstant[
+ cyl_rx = {150 , Name StrCat[pp0 , "1cylinder X-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ cyl_ry = {150 , Name StrCat[pp0 , "2cylinder Y-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ cyl_h = {200 , Name StrCat[pp0 , "3cylinder height [nm]"] , Highlight Str[colorppOK] , Closed 0}];
+ // FIXME gmsh Max?
+ If (cyl_rx>cyl_ry)
+ rbb_tmp=cyl_rx;
+ Else rbb_tmp=cyl_ry;
+ EndIf
+ rbb = 0.5*Sqrt[rbb_tmp^2+cyl_h^2];
+EndIf
+If (flag_shape==CONE)
+ DefineConstant[
+ cone_rx = {300 , Name StrCat[pp0 , "1cone basis X-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ cone_ry = {300 , Name StrCat[pp0 , "1cone basis Y-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ cone_h = {300 , Name StrCat[pp0 , "2cone height [nm]"] , Highlight Str[colorppOK] , Closed 0}];
+ // FIXME gmsh Max?
+ If (cone_rx>cone_ry)
+ rbb_tmp=cone_rx;
+ Else rbb_tmp=cone_ry;
+ EndIf
+ If (2.0*cone_h/3.0>Sqrt[rbb_tmp^2+cone_h^2/9.0])
+ rbb=2.0*cone_h/3.0;
+ Else rbb=Sqrt[rbb_tmp^2+cone_h^2/9.0];
+ EndIf
+EndIf
+If (flag_shape==TOR)
+ DefineConstant[
+ tor_r1 = { 300 , Name StrCat[pp0 , "1torus radius 1 [nm]"], Highlight Str[colorppOK] , Closed 0},
+ tor_r2x = { 100 , Name StrCat[pp0 , "2torus radius 2x [nm]"], Highlight Str[colorppOK] , Closed 0},
+ tor_r2z = { 50 , Name StrCat[pp0 , "3torus radius 2z [nm]"], Highlight Str[colorppOK] , Closed 0},
+ tor_angle = { 340 , Name StrCat[pp0 , "4torus angle [deg]"] , Highlight Str[colorppOK] , Closed 0, Min 5, Max 355}];
+ rbb = tor_r1+tor_r2x;
+EndIf
+DefineConstant[
+ rot_theta = {0 , Name StrCat[pp0 , "5rotate scatterer (polar) [deg]"] , Highlight Str[colorppOK] , Closed 0, Min 0, Max 180},
+ rot_phi = {0 , Name StrCat[pp0 , "6rotate scatterer (azimut) [deg]"], Highlight Str[colorppOK] , Closed 0, Min 0, Max 360}];
+
+flag_study = RES_TMAT;
+
+DefineConstant[
+ epsr_In_re = { 9., Name StrCat[pp2 , "0scatterer permittivity (real) []"] , Highlight Str[colorppOK] , Closed 0},
+ epsr_In_im = { 0., Name StrCat[pp2 , "1scatterer permittivity (imag) []"] , Highlight Str[colorppOK] , Closed 0},
+ epsr_Out_re = { 1., Name StrCat[pp2 , "2background permittivity (real) []"], Highlight Str[colorppOK] , Closed 0},
+ epsr_Out_im = { 0., Name StrCat[pp2 , "3background permittivity (imag) []"], Highlight Str[colorppOK] , Closed 0}
+];
+
+DefineConstant[
+ lambda0 = {587.6 , Name StrCat[pp2 , "4wavelength [nm]"] , Highlight Str[colorppOK] , Closed 0,GmshOption "Reset", Autocheck 0}];
+
+lambda_bg = lambda0/Sqrt[epsr_Out_re];
+
+DefineConstant[
+ n_max = {1 , Name StrCat[pp2 , "8n_max integer"], Highlight Str[colorppOK] , Closed 0, Min 0, Max 5},
+ siwt = {0 , Name StrCat[pp2 , "9Time sign e^(+|-iwt)"], Choices {0="e^(-iwt)",2="e^(+iwt)"} , Closed 0}
+];
+flag_cartpml = 0;
+DefineConstant[
+ space2pml = {lambda_bg/2 , Name StrCat[pp3 , "0space around scatterer [nm]"] , Highlight Str[colorppNO] , Closed 1,Min (lambda_bg/10.), Max (3*lambda_bg)},
+ pml_size = {lambda_bg , Name StrCat[pp3 , "1PML thickness [nm]"] , Highlight Str[colorppWA] , Closed 0, Min (lambda_bg/10), Max (3*lambda_bg)},
+ paramaille = {5. , Name StrCat[pp3 , "2mesh size"], Highlight Str[colorppWA] , Closed 0},
+ refine_scat = {1. , Name StrCat[pp3 , "3scatterer mesh refinement"], Highlight Str[colorppWA] , Closed 0}
+ is_FEM_o2 = {1 , Name StrCat[pp3 , "4Interpolation order "] , Choices {0="order 1",1="order 2"}, Closed 0}
+];
+
+DefineConstant[
+ flag_plotcuts = {1, Choices{0,1}, Name StrCat[pp5, "Plot radial cuts?"]}
+];
+
+
+
+// FIXME conditional for normalization
+If (flag_shape==ELL)
+ ell_rx = ell_rx*nm;
+ ell_ry = ell_ry*nm;
+ ell_rz = ell_rz*nm;
+EndIf
+If (flag_shape==PARALL)
+ par_ax = par_ax*nm;
+ par_ay = par_ay*nm;
+ par_az = par_az*nm;
+EndIf
+If (flag_shape==CYL)
+ cyl_rx = cyl_rx*nm;
+ cyl_ry = cyl_ry*nm;
+ cyl_h = cyl_h *nm;
+EndIf
+If (flag_shape==CONE)
+ cone_rx = cone_rx*nm;
+ cone_ry = cone_ry*nm;
+ cone_h = cone_h*nm;
+EndIf
+If (flag_shape==TOR)
+ tor_r1 = tor_r1*nm;
+ tor_r2x = tor_r2x*nm;
+ tor_r2z = tor_r2z*nm;
+EndIf
+
+lambda0 = lambda0*nm;
+lambda_bg = lambda_bg*nm;
+pml_size = pml_size*nm;
+space2pml = space2pml*nm;
+rbb = rbb*nm;
+
+r_pml_in = space2pml+rbb;
+r_pml_out = space2pml+rbb+pml_size;
+
+
+p_max = n_max*n_max+2*n_max;
+siwt=siwt-1;
+
+sph_scan = 0.00001;
+
+npts_plot_theta = 25;
+npts_plot_phi = 50;