new test case with 4 spheres

ElectromagneticScattering/scattererTmatrix.geo

@@ -31,18 +31,29 @@ If (flag_shape==CONE)
   Dilate { { 0,0,0 }, { 1 , cone_ry/cone_rx , 1} } { Volume{1}; }
 EndIf
 If (flag_shape==TOR)
-  U (1) = {0,0,0,tor_r1,tor_r2x,tor_angle*Pi/180};
+  Torus(1) = {0,0,0,tor_r1,tor_r2x,tor_angle*Pi/180};
   Dilate { { 0,0,0 }, { 1 , 1 , tor_r2z/tor_r2x} } { Volume{1}; }
 EndIf
+If (flag_shape==QUADSPH)
+  Sphere(4) = {x_sph1,y_sph1,z_sph1,r_sph1};
+  Sphere(5) = {x_sph2,y_sph2,z_sph2,r_sph2};
+  Sphere(6) = {x_sph3,y_sph3,z_sph3,r_sph3};
+  Sphere(7) = {x_sph4,y_sph4,z_sph4,r_sph4};  
+EndIf
 
 Sphere (2) = { 0,0,0,r_pml_in}; 
 Sphere (3) = { 0,0,0,r_pml_out}; 
 
 Coherence;
-
+If (flag_shape==QUADSPH)
+  Physical Volume("Scatterer" ,1) = {4,5,6,7};
+  Physical Volume("Background",2) = {8};
+  Physical Volume("PML",       3) = {9};
+Else
   Physical Volume("Scatterer" ,1) = {1};
   Physical Volume("Background",2) = {2};
   Physical Volume("PML"       ,3) = {3};
+EndIf
 
 If(flag_shape==ELL)
   Physical Surface("SurfInt",10) = {2};
@@ -64,6 +75,10 @@ If(flag_shape==TOR)
   Physical Surface("SurfInt",10) = {4};
   Physical Surface("SurfPML",20) = {5};
 EndIf
+If(flag_shape==QUADSPH)
+  Physical Surface("SurfInt",10) = {5};
+  Physical Surface("SurfPML",20) = {6};
+EndIf
 
 Characteristic Length{PointsOf{Physical Volume{3};}} = PML_lc;
 Characteristic Length{PointsOf{Physical Volume{2};}} = Out_lc;

ElectromagneticScattering/scattererTmatrix_data.geo

@@ -26,6 +26,7 @@ PARALL = 1;
 CYL     = 2;
 CONE    = 3;
 TOR     = 4;
+QUADSPH = 5;
 
 RES_PW    = 0;
 RES_TMAT  = 1;
@@ -34,7 +35,7 @@ RES_QNM   = 3;
 
 DefineConstant[
   flag_shape = {ELL , Name StrCat[pp0, "0Scatterer shape"], 
-    Choices {ELL="ellispoid",PARALL="parallelepiped",CYL="cylinder",CONE="cone",TOR="split U"},Closed 0}];
+    Choices {ELL="ellispoid",PARALL="parallelepiped",CYL="cylinder",CONE="cone",TOR="split U", QUADSPH="4 spheres"},Closed 0}];
 
 If (flag_shape==ELL)
   DefineConstant[
@@ -94,6 +95,16 @@ If (flag_shape==TOR)
     tor_angle = { 340 , Name StrCat[pp0 , "4U angle [deg]"]   , Highlight Str[colorppOK] , Closed 0, Min 5, Max 355}];
   rbb = tor_r1+tor_r2x;
 EndIf
+If (flag_shape==QUADSPH)
+  DefineConstant[
+    r_sph1 = {50.0 , Name "sphere 1 radius [nm]" },
+    r_sph2 = {60.0 , Name "sphere 2 radius [nm]" },
+    r_sph3 = {70.0 , Name "sphere 3 radius [nm]" },
+    r_sph4 = {80.0 , Name "sphere 4 radius [nm]" },
+    a = {300.0 , Name "tetrahedron side length [nm]" }
+  ];
+  rbb = 0.6*a;
+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}];
@@ -157,6 +168,25 @@ If (flag_shape==TOR)
   tor_r2x = tor_r2x*nm;
   tor_r2z = tor_r2z*nm;
 EndIf
+If (flag_shape==QUADSPH)
+  a      = a*nm;
+  r_sph1 = r_sph1*nm;
+  r_sph2 = r_sph2*nm;
+  r_sph3 = r_sph3*nm;
+  r_sph4 = r_sph4*nm;
+  x_sph1 = -a/2;
+  y_sph1 = -a*Sqrt[3]/6;
+  z_sph1 = -a*Sqrt[6]/12;
+  x_sph2 = a/2;
+  y_sph2 = -a*Sqrt[3]/6;
+  z_sph2 = -a*Sqrt[6]/12;
+  x_sph3 = 0;
+  y_sph3 = a*Sqrt[3]/3;
+  z_sph3 = -a*Sqrt[6]/12;
+  x_sph4 = 0;
+  y_sph4 = 0;
+  z_sph4 = a*Sqrt[6]/4;
+EndIf
 
 lambda0   = lambda0*nm;
 lambda_bg = lambda_bg*nm;