diff --git a/ElectromagneticScattering/README.md b/ElectromagneticScattering/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..2fd79dd011cffae33b06e2db06ebfc1a24fd81ab
--- /dev/null
+++ b/ElectromagneticScattering/README.md
@@ -0,0 +1,23 @@
+A Onelab model for 3D scattering problems in nanophotonics.
+
+## Synopsis
+
+This project contains a [Onelab](http://onelab.info/wiki/ONELAB) model for solving general 3D electromagnetic scattering problems:
+* T-matrix computation: See https://arxiv.org/abs/1802.00596 for details
+* Quasi-normal modes
+* Plane wave response
+* Green's function (LDOS)
+
+## Installation
+
+This model requires the following (open-source) programs:
+* [onelab](http://onelab.info/wiki/ONELAB), which bundles both [gmsh](http://www.gmsh.info/) and [getdp](http://www.getdp.info/)
+* python (v2.7.x or v3.6.x)
+
+## Running the model
+
+Open `scattering.pro` with Gmsh.
+
+## Authors
+
+Guillaume Demésy
\ No newline at end of file
diff --git a/ElectromagneticScattering/scattering.geo b/ElectromagneticScattering/scattering.geo
new file mode 100644
index 0000000000000000000000000000000000000000..b51d28eb4bfe11e95a8a4c9d1f13cc6d1fbd812c
--- /dev/null
+++ b/ElectromagneticScattering/scattering.geo
@@ -0,0 +1,753 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering.geo //
+///////////////////////////////
+
+Include "scattering_data.geo";
+In_n = Sqrt[Fabs[epsr_In_re]];
+// Out_n = Sqrt[Fabs[epsr_Out_re]];
+
+paramaille_pml = paramaille/1.1;
+
+Out_lc = lambda_bg/paramaille;
+PML_lc = lambda_bg/paramaille_pml;
+In_lc = lambda0/(paramaille*In_n*refine_scat);
+CenterScat_lc = lambda0/(paramaille*In_n);
+
+If(flag_cartpml==1)
+ dom_x = r_pml_in*2;
+ dom_y = r_pml_in*2;
+ dom_z = r_pml_in*2;
+ PML_bot = pml_size;
+ PML_lat = pml_size;
+ PML_top = pml_size;
+ Point(1) = {-dom_x/2,-dom_y/2, -dom_z/2, Out_lc};
+ Point(2) = {-dom_x/2, dom_y/2, -dom_z/2, Out_lc};
+ Point(3) = { dom_x/2, dom_y/2, -dom_z/2, Out_lc};
+ Point(4) = { dom_x/2,-dom_y/2, -dom_z/2, Out_lc};
+ Point(5) = {-dom_x/2,-dom_y/2, dom_z/2, Out_lc};
+ Point(6) = {-dom_x/2, dom_y/2, dom_z/2, Out_lc};
+ Point(7) = { dom_x/2, dom_y/2, dom_z/2, Out_lc};
+ Point(8) = { dom_x/2,-dom_y/2, dom_z/2, Out_lc};
+ Point(9) = {-dom_x/2,-dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(10) = {-dom_x/2, dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(11) = { dom_x/2, dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(12) = { dom_x/2,-dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(13) = {-dom_x/2,-dom_y/2, dom_z/2+PML_top, PML_lc};
+ Point(14) = {-dom_x/2, dom_y/2, dom_z/2+PML_top, PML_lc};
+ Point(15) = { dom_x/2, dom_y/2, dom_z/2+PML_top, PML_lc};
+ Point(16) = { dom_x/2,-dom_y/2, dom_z/2+PML_top, PML_lc};
+ Point(17) = {-dom_x/2,-dom_y/2-PML_lat, -dom_z/2 , PML_lc};
+ Point(18) = {-dom_x/2, dom_y/2+PML_lat, -dom_z/2 , PML_lc};
+ Point(19) = { dom_x/2, dom_y/2+PML_lat, -dom_z/2 , PML_lc};
+ Point(20) = { dom_x/2,-dom_y/2-PML_lat, -dom_z/2 , PML_lc};
+ Point(21) = {-dom_x/2,-dom_y/2-PML_lat, dom_z/2 , PML_lc};
+ Point(22) = {-dom_x/2, dom_y/2+PML_lat, dom_z/2 , PML_lc};
+ Point(23) = { dom_x/2, dom_y/2+PML_lat, dom_z/2 , PML_lc};
+ Point(24) = { dom_x/2,-dom_y/2-PML_lat, dom_z/2 , PML_lc};
+ Point(25) = {-dom_x/2-PML_lat,-dom_y/2, -dom_z/2 , PML_lc};
+ Point(26) = {-dom_x/2-PML_lat, dom_y/2, -dom_z/2 , PML_lc};
+ Point(27) = { dom_x/2+PML_lat, dom_y/2, -dom_z/2 , PML_lc};
+ Point(28) = { dom_x/2+PML_lat,-dom_y/2, -dom_z/2 , PML_lc};
+ Point(29) = {-dom_x/2-PML_lat,-dom_y/2, dom_z/2 , PML_lc};
+ Point(30) = {-dom_x/2-PML_lat, dom_y/2, dom_z/2 , PML_lc};
+ Point(31) = { dom_x/2+PML_lat, dom_y/2, dom_z/2 , PML_lc};
+ Point(32) = { dom_x/2+PML_lat,-dom_y/2, dom_z/2 , PML_lc};
+ Point(33) = {-dom_x/2-PML_lat,-dom_y/2-PML_lat, -dom_z/2-PML_bot, PML_lc};
+ Point(34) = {-dom_x/2-PML_lat, dom_y/2+PML_lat, -dom_z/2-PML_bot, PML_lc};
+ Point(35) = { dom_x/2+PML_lat, dom_y/2+PML_lat, -dom_z/2-PML_bot, PML_lc};
+ Point(36) = { dom_x/2+PML_lat,-dom_y/2-PML_lat, -dom_z/2-PML_bot, PML_lc};
+ Point(37) = {-dom_x/2-PML_lat,-dom_y/2-PML_lat, dom_z/2+PML_top, PML_lc};
+ Point(38) = {-dom_x/2-PML_lat, dom_y/2+PML_lat, dom_z/2+PML_top, PML_lc};
+ Point(39) = { dom_x/2+PML_lat, dom_y/2+PML_lat, dom_z/2+PML_top, PML_lc};
+ Point(40) = { dom_x/2+PML_lat,-dom_y/2-PML_lat, dom_z/2+PML_top, PML_lc};
+ Point(41) = {-dom_x/2,-dom_y/2-PML_lat, -dom_z/2-PML_bot, PML_lc};
+ Point(42) = {-dom_x/2, dom_y/2+PML_lat, -dom_z/2-PML_bot, PML_lc};
+ Point(43) = { dom_x/2,-dom_y/2-PML_lat, -dom_z/2-PML_top, PML_lc};
+ Point(44) = { dom_x/2, dom_y/2+PML_lat, -dom_z/2-PML_top, PML_lc};
+ Point(45) = {-dom_x/2,-dom_y/2-PML_lat, dom_z/2+PML_bot, PML_lc};
+ Point(46) = {-dom_x/2, dom_y/2+PML_lat, dom_z/2+PML_bot, PML_lc};
+ Point(47) = { dom_x/2,-dom_y/2-PML_lat, dom_z/2+PML_top, PML_lc};
+ Point(48) = { dom_x/2, dom_y/2+PML_lat, dom_z/2+PML_top, PML_lc};
+ Point(49) = {-dom_x/2-PML_lat, -dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(50) = { dom_x/2+PML_lat, -dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(51) = {-dom_x/2-PML_lat, -dom_y/2, dom_z/2+PML_bot, PML_lc};
+ Point(52) = { dom_x/2+PML_lat, -dom_y/2, dom_z/2+PML_bot, PML_lc};
+ Point(53) = {-dom_x/2-PML_lat, dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(54) = { dom_x/2+PML_lat, dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(55) = {-dom_x/2-PML_lat, dom_y/2, dom_z/2+PML_bot, PML_lc};
+ Point(56) = { dom_x/2+PML_lat, dom_y/2, dom_z/2+PML_bot, PML_lc};
+ Point(57) = {-dom_x/2-PML_lat, -dom_y/2-PML_lat, -dom_z/2, PML_lc};
+ Point(58) = { dom_x/2+PML_lat, -dom_y/2-PML_lat, -dom_z/2, PML_lc};
+ Point(59) = {-dom_x/2-PML_lat, dom_y/2+PML_lat, -dom_z/2, PML_lc};
+ Point(60) = { dom_x/2+PML_lat, dom_y/2+PML_lat, -dom_z/2, PML_lc};
+ Point(61) = {-dom_x/2-PML_lat, -dom_y/2-PML_lat, dom_z/2, PML_lc};
+ Point(62) = { dom_x/2+PML_lat, -dom_y/2-PML_lat, dom_z/2, PML_lc};
+ Point(63) = {-dom_x/2-PML_lat, dom_y/2+PML_lat, dom_z/2, PML_lc};
+ Point(64) = { dom_x/2+PML_lat, dom_y/2+PML_lat, dom_z/2, PML_lc};
+ Point(73) = { 0,0,0, CenterScat_lc};
+
+ Line(1) = {33, 49};
+ Line(2) = {49, 53};
+ Line(3) = {53, 34};
+ Line(4) = {34, 42};
+ Line(5) = {42, 44};
+ Line(6) = {44, 35};
+ Line(7) = {35, 54};
+ Line(8) = {54, 11};
+ Line(9) = {11, 10};
+ Line(10) = {10, 53};
+ Line(11) = {49, 9};
+ Line(12) = {9, 12};
+ Line(13) = {12, 50};
+ Line(14) = {50, 54};
+ Line(15) = {50, 36};
+ Line(16) = {36, 43};
+ Line(17) = {43, 41};
+ Line(18) = {41, 33};
+ Line(19) = {41, 9};
+ Line(20) = {9, 10};
+ Line(21) = {10, 42};
+ Line(22) = {43, 12};
+ Line(23) = {12, 11};
+ Line(24) = {11, 44};
+ Line(25) = {57, 25};
+ Line(26) = {25, 26};
+ Line(27) = {26, 59};
+ Line(28) = {59, 18};
+ Line(29) = {18, 19};
+ Line(30) = {19, 60};
+ Line(31) = {60, 27};
+ Line(32) = {27, 28};
+ Line(33) = {28, 58};
+ Line(34) = {58, 20};
+ Line(35) = {20, 17};
+ Line(36) = {17, 57};
+ Line(37) = {25, 1};
+ Line(38) = {1, 4};
+ Line(39) = {4, 28};
+ Line(40) = {27, 3};
+ Line(41) = {3, 2};
+ Line(42) = {2, 26};
+ Line(43) = {17, 1};
+ Line(44) = {1, 2};
+ Line(45) = {2, 18};
+ Line(46) = {19, 3};
+ Line(47) = {3, 4};
+ Line(48) = {4, 20};
+ Line(49) = {61, 21};
+ Line(50) = {21, 24};
+ Line(51) = {24, 62};
+ Line(52) = {62, 32};
+ Line(53) = {32, 31};
+ Line(54) = {31, 64};
+ Line(55) = {64, 23};
+ Line(56) = {23, 22};
+ Line(57) = {22, 63};
+ Line(58) = {63, 30};
+ Line(59) = {30, 29};
+ Line(60) = {29, 61};
+ Line(61) = {21, 5};
+ Line(62) = {5, 6};
+ Line(63) = {6, 22};
+ Line(64) = {23, 7};
+ Line(65) = {7, 8};
+ Line(66) = {8, 24};
+ Line(67) = {32, 8};
+ Line(68) = {8, 5};
+ Line(69) = {5, 29};
+ Line(70) = {30, 6};
+ Line(71) = {6, 7};
+ Line(72) = {7, 31};
+ Line(73) = {37, 45};
+ Line(74) = {45, 47};
+ Line(75) = {47, 40};
+ Line(76) = {40, 52};
+ Line(77) = {52, 56};
+ Line(78) = {56, 39};
+ Line(79) = {39, 48};
+ Line(80) = {48, 46};
+ Line(81) = {46, 38};
+ Line(82) = {38, 55};
+ Line(83) = {55, 51};
+ Line(84) = {51, 37};
+ Line(85) = {45, 13};
+ Line(86) = {13, 14};
+ Line(87) = {14, 46};
+ Line(88) = {55, 14};
+ Line(89) = {14, 15};
+ Line(90) = {15, 56};
+ Line(91) = {52, 16};
+ Line(92) = {16, 13};
+ Line(93) = {13, 51};
+ Line(94) = {48, 15};
+ Line(95) = {15, 16};
+ Line(96) = {16, 47};
+ Line(97) = {37, 61};
+ Line(98) = {61, 57};
+ Line(99) = {57, 33};
+ Line(100) = {45, 21};
+ Line(101) = {21, 17};
+ Line(102) = {17, 41};
+ Line(103) = {47, 24};
+ Line(104) = {24, 20};
+ Line(105) = {20, 43};
+ Line(106) = {40, 62};
+ Line(107) = {62, 58};
+ Line(108) = {58, 36};
+ Line(109) = {52, 32};
+ Line(110) = {32, 28};
+ Line(111) = {28, 50};
+ Line(112) = {56, 31};
+ Line(113) = {31, 27};
+ Line(114) = {27, 54};
+ Line(115) = {39, 64};
+ Line(116) = {64, 60};
+ Line(117) = {60, 35};
+ Line(118) = {48, 23};
+ Line(119) = {23, 19};
+ Line(120) = {19, 44};
+ Line(121) = {46, 22};
+ Line(122) = {22, 18};
+ Line(123) = {18, 42};
+ Line(124) = {38, 63};
+ Line(125) = {63, 59};
+ Line(126) = {59, 34};
+ Line(127) = {55, 30};
+ Line(128) = {30, 26};
+ Line(129) = {26, 53};
+ Line(130) = {51, 29};
+ Line(131) = {29, 25};
+ Line(132) = {25, 49};
+ Line(133) = {14, 6};
+ Line(134) = {6, 2};
+ Line(135) = {2, 10};
+ Line(136) = {15, 7};
+ Line(137) = {7, 3};
+ Line(138) = {3, 11};
+ Line(139) = {16, 8};
+ Line(140) = {8, 4};
+ Line(141) = {4, 12};
+ Line(142) = {13, 5};
+ Line(143) = {5, 1};
+ Line(144) = {1, 9};
+
+ Line Loop(173) = {80, -87, 89, -94}; Plane Surface(174) = {173};
+ Line Loop(175) = {94, 90, 78, 79}; Plane Surface(176) = {175};
+ Line Loop(177) = {77, -90, 95, -91}; Plane Surface(178) = {177};
+ Line Loop(179) = {95, 92, 86, 89}; Plane Surface(180) = {179};
+ Line Loop(181) = {91, 96, 75, 76}; Plane Surface(182) = {181};
+ Line Loop(183) = {92, -85, 74, -96}; Plane Surface(184) = {183};
+ Line Loop(185) = {82, 88, 87, 81}; Plane Surface(186) = {185};
+ Line Loop(187) = {86, -88, 83, -93}; Plane Surface(188) = {187};
+ Line Loop(189) = {93, 84, 73, 85}; Plane Surface(190) = {189};
+ Line Loop(191) = {58, 70, 63, 57}; Plane Surface(192) = {191};
+ Line Loop(193) = {63, -56, 64, -71}; Plane Surface(194) = {193};
+ Line Loop(195) = {64, 72, 54, 55}; Plane Surface(196) = {195};
+ Line Loop(197) = {72, -53, 67, -65}; Plane Surface(198) = {197};
+ Line Loop(199) = {67, 66, 51, 52}; Plane Surface(200) = {199};
+ Line Loop(201) = {71, 65, 68, 62}; Plane Surface(202) = {201};
+ Line Loop(203) = {68, -61, 50, -66}; Plane Surface(204) = {203};
+ Line Loop(205) = {59, -69, 62, -70}; Plane Surface(206) = {205};
+ Line Loop(207) = {69, 60, 49, 61}; Plane Surface(208) = {207};
+ Line Loop(209) = {28, -45, 42, 27}; Plane Surface(210) = {209};
+ Line Loop(211) = {29, 46, 41, 45}; Plane Surface(212) = {211};
+ Line Loop(213) = {30, 31, 40, -46}; Plane Surface(214) = {213};
+ Line Loop(215) = {40, 47, 39, -32}; Plane Surface(216) = {215};
+ Line Loop(217) = {41, -44, 38, -47}; Plane Surface(218) = {217};
+ Line Loop(219) = {39, 33, 34, -48}; Plane Surface(220) = {219};
+ Line Loop(221) = {38, 48, 35, 43}; Plane Surface(222) = {221};
+ Line Loop(223) = {26, -42, -44, -37}; Plane Surface(224) = {223};
+ Line Loop(225) = {37, -43, 36, 25}; Plane Surface(226) = {225};
+ Line Loop(227) = {3, 4, -21, 10}; Plane Surface(228) = {227};
+ Line Loop(229) = {5, -24, 9, 21}; Plane Surface(230) = {229};
+ Line Loop(231) = {6, 7, 8, 24}; Plane Surface(232) = {231};
+ Line Loop(233) = {14, 8, -23, 13}; Plane Surface(234) = {233};
+ Line Loop(235) = {9, -20, 12, 23}; Plane Surface(236) = {235};
+ Line Loop(237) = {2, -10, -20, -11}; Plane Surface(238) = {237};
+ Line Loop(239) = {13, 15, 16, 22}; Plane Surface(240) = {239};
+ Line Loop(241) = {12, -22, 17, 19}; Plane Surface(242) = {241};
+ Line Loop(243) = {19, -11, -1, -18}; Plane Surface(244) = {243};
+ Line Loop(245) = {127, 70, -133, -88}; Plane Surface(246) = {245};
+ Line Loop(247) = {133, 71, -136, -89}; Plane Surface(248) = {247};
+ Line Loop(249) = {90, 112, -72, -136}; Plane Surface(250) = {249};
+ Line Loop(251) = {128, -42, -134, -70}; Plane Surface(252) = {251};
+ Line Loop(253) = {41, -134, 71, 137}; Plane Surface(254) = {253};
+ Line Loop(255) = {137, -40, -113, -72}; Plane Surface(256) = {255};
+ Line Loop(257) = {129, -10, -135, 42}; Plane Surface(258) = {257};
+ Line Loop(259) = {135, -9, -138, 41}; Plane Surface(260) = {259};
+ Line Loop(261) = {138, -8, -114, 40}; Plane Surface(262) = {261};
+ Line Loop(263) = {130, -69, -142, 93}; Plane Surface(264) = {263};
+ Line Loop(265) = {92, 142, -68, -139}; Plane Surface(266) = {265};
+ Line Loop(267) = {139, -67, -109, 91}; Plane Surface(268) = {267};
+ Line Loop(269) = {131, 37, -143, 69}; Plane Surface(270) = {269};
+ Line Loop(271) = {143, 38, -140, 68}; Plane Surface(272) = {271};
+ Line Loop(273) = {140, 39, -110, 67}; Plane Surface(274) = {273};
+ Line Loop(275) = {132, 11, -144, -37}; Plane Surface(276) = {275};
+ Line Loop(277) = {12, -141, -38, 144}; Plane Surface(278) = {277};
+ Line Loop(279) = {141, 13, -111, -39}; Plane Surface(280) = {279};
+ Line Loop(281) = {99, -18, -102, 36}; Plane Surface(282) = {281};
+ Line Loop(283) = {102, -17, -105, 35}; Plane Surface(284) = {283};
+ Line Loop(285) = {34, 105, -16, -108}; Plane Surface(286) = {285};
+ Line Loop(287) = {34, -104, 51, 107}; Plane Surface(288) = {287};
+ Line Loop(289) = {50, 104, 35, -101}; Plane Surface(290) = {289};
+ Line Loop(291) = {36, -98, 49, 101}; Plane Surface(292) = {291};
+ Line Loop(293) = {97, 49, -100, -73}; Plane Surface(294) = {293};
+ Line Loop(295) = {74, 103, -50, -100}; Plane Surface(296) = {295};
+ Line Loop(297) = {103, 51, -106, -75}; Plane Surface(298) = {297};
+ Line Loop(299) = {126, 4, -123, -28}; Plane Surface(300) = {299};
+ Line Loop(301) = {125, 28, -122, 57}; Plane Surface(302) = {301};
+ Line Loop(303) = {122, 29, -119, 56}; Plane Surface(304) = {303};
+ Line Loop(305) = {119, 30, -116, 55}; Plane Surface(306) = {305};
+ Line Loop(307) = {55, -118, -79, 115}; Plane Surface(308) = {307};
+ Line Loop(309) = {80, 121, -56, -118}; Plane Surface(310) = {309};
+ Line Loop(311) = {81, 124, -57, -121}; Plane Surface(312) = {311};
+ Line Loop(313) = {123, 5, -120, -29}; Plane Surface(314) = {313};
+ Line Loop(315) = {120, 6, -117, -30}; Plane Surface(316) = {315};
+ Line Loop(317) = {121, -63, -133, 87}; Plane Surface(318) = {317};
+ Line Loop(319) = {133, -62, -142, 86}; Plane Surface(320) = {319};
+ Line Loop(321) = {142, -61, -100, 85}; Plane Surface(322) = {321};
+ Line Loop(323) = {122, -45, -134, 63}; Plane Surface(324) = {323};
+ Line Loop(325) = {134, -44, -143, 62}; Plane Surface(326) = {325};
+ Line Loop(327) = {143, -43, -101, 61}; Plane Surface(328) = {327};
+ Line Loop(329) = {123, -21, -135, 45}; Plane Surface(330) = {329};
+ Line Loop(331) = {135, -20, -144, 44}; Plane Surface(332) = {331};
+ Line Loop(333) = {144, -19, -102, 43}; Plane Surface(334) = {333};
+ Line Loop(335) = {126, -3, -129, 27}; Plane Surface(336) = {335};
+ Line Loop(337) = {129, -2, -132, 26}; Plane Surface(338) = {337};
+ Line Loop(339) = {132, -1, -99, 25}; Plane Surface(340) = {339};
+ Line Loop(341) = {125, -27, -128, -58}; Plane Surface(342) = {341};
+ Line Loop(343) = {128, -26, -131, -59}; Plane Surface(344) = {343};
+ Line Loop(345) = {131, -25, -98, -60}; Plane Surface(346) = {345};
+ Line Loop(347) = {97, -60, -130, 84}; Plane Surface(348) = {347};
+ Line Loop(349) = {83, 130, -59, -127}; Plane Surface(350) = {349};
+ Line Loop(351) = {58, -127, -82, 124}; Plane Surface(352) = {351};
+ Line Loop(353) = {96, 103, -66, -139}; Plane Surface(354) = {353};
+ Line Loop(355) = {139, -65, -136, 95}; Plane Surface(356) = {355};
+ Line Loop(357) = {136, -64, -118, 94}; Plane Surface(358) = {357};
+ Line Loop(359) = {66, 104, -48, -140}; Plane Surface(360) = {359};
+ Line Loop(361) = {140, -47, -137, 65}; Plane Surface(362) = {361};
+ Line Loop(363) = {137, -46, -119, 64}; Plane Surface(364) = {363};
+ Line Loop(365) = {105, 22, -141, 48}; Plane Surface(366) = {365};
+ Line Loop(367) = {141, 23, -138, 47}; Plane Surface(368) = {367};
+ Line Loop(369) = {138, 24, -120, 46}; Plane Surface(370) = {369};
+ Line Loop(371) = {108, -15, -111, 33}; Plane Surface(372) = {371};
+ Line Loop(373) = {111, 14, -114, 32}; Plane Surface(374) = {373};
+ Line Loop(375) = {114, -7, -117, 31}; Plane Surface(376) = {375};
+ Line Loop(377) = {107, -33, -110, -52}; Plane Surface(378) = {377};
+ Line Loop(379) = {110, -32, -113, -53}; Plane Surface(380) = {379};
+ Line Loop(381) = {113, -31, -116, -54}; Plane Surface(382) = {381};
+ Line Loop(383) = {106, 52, -109, -76}; Plane Surface(384) = {383};
+ Line Loop(385) = {109, 53, -112, -77}; Plane Surface(386) = {385};
+ Line Loop(387) = {112, 54, -115, -78}; Plane Surface(388) = {387};
+
+ Surface Loop(393) = {182, 298, 384, 354, 200, 268}; Volume(394) = {393};
+ Surface Loop(395) = {184, 296, 354, 266, 204, 322}; Volume(396) = {395};
+ Surface Loop(397) = {190, 348, 294, 264, 208, 322}; Volume(398) = {397};
+ Surface Loop(399) = {188, 350, 264, 246, 320, 206}; Volume(400) = {399};
+ Surface Loop(401) = {186, 352, 312, 192, 318, 246}; Volume(402) = {401};
+ Surface Loop(403) = {180, 248, 266, 202, 356, 320}; Volume(404) = {403};
+ Surface Loop(405) = {310, 174, 318, 248, 194, 358}; Volume(406) = {405};
+ Surface Loop(407) = {176, 388, 308, 358, 250, 196}; Volume(408) = {407};
+ Surface Loop(409) = {178, 386, 198, 356, 268, 250}; Volume(410) = {409};
+ Surface Loop(411) = {378, 288, 200, 274, 360, 220}; Volume(412) = {411};
+ Surface Loop(413) = {216, 380, 198, 274, 256, 362}; Volume(414) = {413};
+ Surface Loop(415) = {306, 382, 214, 364, 196, 256}; Volume(416) = {415};
+ Surface Loop(417) = {304, 194, 364, 212, 254, 324}; Volume(418) = {417};
+ Surface Loop(419) = {302, 342, 192, 324, 252, 210}; Volume(420) = {419};
+ Surface Loop(421) = {344, 206, 252, 224, 270, 326}; Volume(422) = {421};
+ Surface Loop(423) = {292, 346, 328, 208, 270, 226}; Volume(424) = {423};
+ Surface Loop(425) = {290, 360, 204, 272, 328, 222}; Volume(426) = {425};
+ Surface Loop(427) = {240, 372, 286, 220, 366, 280}; Volume(428) = {427};
+ Surface Loop(429) = {234, 374, 216, 368, 262, 280}; Volume(430) = {429};
+ Surface Loop(431) = {376, 232, 316, 370, 214, 262}; Volume(432) = {431};
+ Surface Loop(433) = {314, 230, 260, 212, 330, 370}; Volume(434) = {433};
+ Surface Loop(435) = {228, 336, 300, 210, 258, 330}; Volume(436) = {435};
+ Surface Loop(437) = {236, 332, 368, 278, 260, 218}; Volume(438) = {437};
+ Surface Loop(439) = {284, 242, 366, 334, 222, 278}; Volume(440) = {439};
+ Surface Loop(441) = {282, 340, 244, 334, 226, 276}; Volume(442) = {441};
+ Surface Loop(443) = {238, 338, 258, 276, 332, 224}; Volume(444) = {443};
+
+ pt_sc=newp;
+ // Printf("newp scatterer=%g",pt_sc);
+ l_sc =newl;
+ ll_sc =newll;
+ s_sc =news;
+
+ If(flag_shape==ELL)
+ Point(pt_sc+20) = { 0, 0, ell_rz, In_lc}; Point(pt_sc+21) = { 0, 0,-ell_rz, In_lc};
+ Point(pt_sc+22) = { 0, ell_ry, 0, In_lc}; Point(pt_sc+23) = { 0,-ell_ry, 0, In_lc};
+ Point(pt_sc+25) = { ell_rx, 0, 0, In_lc}; Point(pt_sc+26) = {-ell_rx, 0, 0, In_lc};
+ Ellipse(l_sc+145) = {pt_sc+22, 73, pt_sc+26, pt_sc+26} Plane{0,0,1};
+ Ellipse(l_sc+146) = {pt_sc+26, 73, pt_sc+23, pt_sc+23} Plane{0,0,1};
+ Ellipse(l_sc+147) = {pt_sc+23, 73, pt_sc+25, pt_sc+25} Plane{0,0,1};
+ Ellipse(l_sc+148) = {pt_sc+25, 73, pt_sc+22, pt_sc+22} Plane{0,0,1};
+ Ellipse(l_sc+149) = {pt_sc+21, 73, pt_sc+26, pt_sc+26} Plane{0,0,1};
+ Ellipse(l_sc+150) = {pt_sc+26, 73, pt_sc+20, pt_sc+20} Plane{0,0,1};
+ Ellipse(l_sc+151) = {pt_sc+20, 73, pt_sc+25, pt_sc+25} Plane{0,0,1};
+ Ellipse(l_sc+152) = {pt_sc+25, 73, pt_sc+21, pt_sc+21} Plane{0,0,1};
+ Ellipse(l_sc+153) = {pt_sc+23, 73, pt_sc+20, pt_sc+20} Plane{0,0,1};
+ Ellipse(l_sc+154) = {pt_sc+20, 73, pt_sc+22, pt_sc+22} Plane{0,0,1};
+ Ellipse(l_sc+155) = {pt_sc+22, 73, pt_sc+21, pt_sc+21} Plane{0,0,1};
+ Ellipse(l_sc+156) = {pt_sc+21, 73, pt_sc+23, pt_sc+23} Plane{0,0,1};
+ Line Loop(ll_sc+157) = {l_sc+151, -147-l_sc, 153+l_sc}; Surface(158+s_sc) = {157+ll_sc};
+ Line Loop(ll_sc+159) = {l_sc+151, 148+l_sc, -154-l_sc}; Surface(160+s_sc) = {159+ll_sc};
+ Line Loop(ll_sc+161) = {l_sc+145, 150+l_sc, 154+l_sc}; Surface(162+s_sc) = {161+ll_sc};
+ Line Loop(ll_sc+163) = {l_sc+150, -153-l_sc, -146-l_sc}; Surface(164+s_sc) = {163+ll_sc};
+ Line Loop(ll_sc+165) = {l_sc+156, 147+l_sc, 152+l_sc}; Surface(166+s_sc) = {165+ll_sc};
+ Line Loop(ll_sc+167) = {l_sc+152, -155-l_sc, -148-l_sc}; Surface(168+s_sc) = {167+ll_sc};
+ Line Loop(ll_sc+169) = {l_sc+155, 149+l_sc, -145-l_sc}; Surface(170+s_sc) = {169+ll_sc};
+ Line Loop(ll_sc+171) = {l_sc+149, 146+l_sc, -156-l_sc}; Surface(172+s_sc) = {171+ll_sc};
+ Surface Loop(1072) = {160+s_sc, 162+s_sc, 164+s_sc, 158+s_sc, 166+s_sc, 170+s_sc, 172+s_sc, 168+s_sc};
+ Volume(1092) = {1072};
+ EndIf
+ If(flag_shape==PARALL)
+ Point(pt_sc+20) = {-par_ax/2,-par_ay/2,-par_az/2,In_lc};
+ Point(pt_sc+21) = { par_ax/2,-par_ay/2,-par_az/2,In_lc};
+ Point(pt_sc+22) = { par_ax/2, par_ay/2,-par_az/2,In_lc};
+ Point(pt_sc+23) = {-par_ax/2, par_ay/2,-par_az/2,In_lc};
+ Point(pt_sc+24) = {-par_ax/2,-par_ay/2, par_az/2,In_lc};
+ Point(pt_sc+25) = { par_ax/2,-par_ay/2, par_az/2,In_lc};
+ Point(pt_sc+26) = { par_ax/2, par_ay/2, par_az/2,In_lc};
+ Point(pt_sc+27) = {-par_ax/2, par_ay/2, par_az/2,In_lc};
+ Line(l_sc+1) = {pt_sc+20, pt_sc+24};
+ Line(l_sc+2) = {pt_sc+24, pt_sc+25};
+ Line(l_sc+3) = {pt_sc+25, pt_sc+21};
+ Line(l_sc+4) = {pt_sc+21, pt_sc+20};
+ Line(l_sc+5) = {pt_sc+23, pt_sc+27};
+ Line(l_sc+6) = {pt_sc+27, pt_sc+26};
+ Line(l_sc+7) = {pt_sc+26, pt_sc+22};
+ Line(l_sc+8) = {pt_sc+22, pt_sc+23};
+ Line(l_sc+9) = {pt_sc+23, pt_sc+20};
+ Line(l_sc+10) = {pt_sc+22, pt_sc+21};
+ Line(l_sc+11) = {pt_sc+26, pt_sc+25};
+ Line(l_sc+12) = {pt_sc+27, pt_sc+24};
+ Line Loop(ll_sc+13) = {l_sc+2, -11-l_sc, -6-l_sc, l_sc+12};Plane Surface(s_sc+14) = {ll_sc+13};
+ Line Loop(ll_sc+15) = {l_sc+1, -12-l_sc, -5-l_sc, l_sc+9};Plane Surface(s_sc+16) = {ll_sc+15};
+ Line Loop(ll_sc+17) = {l_sc+3, -10-l_sc, -7-l_sc, l_sc+11};Plane Surface(s_sc+18) = {ll_sc+17};
+ Line Loop(ll_sc+19) = {l_sc+8, 5+l_sc, 6+l_sc, l_sc+7};Plane Surface(s_sc+20) = {ll_sc+19};
+ Line Loop(ll_sc+21) = {l_sc+4, -9-l_sc, -8-l_sc, l_sc+10};Plane Surface(s_sc+22) = {ll_sc+21};
+ Line Loop(ll_sc+23) = {l_sc+1, 2+l_sc, 3+l_sc, l_sc+4};Plane Surface(s_sc+24) = {ll_sc+23};
+ Surface Loop(1072) = {s_sc+16, 24+s_sc, 14+s_sc, 18+s_sc, 22+s_sc, 20+s_sc};
+ Volume(1092) = {1072};
+ EndIf
+ If(flag_shape==CYL)
+ Point(pt_sc+20) = { 0, cyl_ry, -cyl_h/2, In_lc};
+ Point(pt_sc+21) = { 0,-cyl_ry, -cyl_h/2, In_lc};
+ Point(pt_sc+22) = { cyl_rx, 0, -cyl_h/2, In_lc};
+ Point(pt_sc+23) = {-cyl_rx, 0, -cyl_h/2, In_lc};
+ Point(pt_sc+24) = {0, 0, -cyl_h/2, In_lc};
+ Point(pt_sc+25) = { 0, cyl_ry, cyl_h/2, In_lc};
+ Point(pt_sc+26) = { 0,-cyl_ry, cyl_h/2, In_lc};
+ Point(pt_sc+27) = { cyl_rx, 0, cyl_h/2, In_lc};
+ Point(pt_sc+28) = {-cyl_rx, 0, cyl_h/2, In_lc};
+ Point(pt_sc+29) = {0, 0, cyl_h/2, In_lc};
+ Ellipse(l_sc+1) = {pt_sc+23, pt_sc+24, pt_sc+21, pt_sc+21};
+ Ellipse(l_sc+2) = {pt_sc+21, pt_sc+24, pt_sc+22, pt_sc+22};
+ Ellipse(l_sc+3) = {pt_sc+22, pt_sc+24, pt_sc+20, pt_sc+20};
+ Ellipse(l_sc+4) = {pt_sc+20, pt_sc+24, pt_sc+24, pt_sc+23};
+ Ellipse(l_sc+5) = {pt_sc+28, pt_sc+29, pt_sc+26, pt_sc+26};
+ Ellipse(l_sc+6) = {pt_sc+26, pt_sc+29, pt_sc+27, pt_sc+27};
+ Ellipse(l_sc+7) = {pt_sc+27, pt_sc+29, pt_sc+25, pt_sc+25};
+ Ellipse(l_sc+8) = {pt_sc+25, pt_sc+29, pt_sc+28, pt_sc+28};
+ Line(l_sc+9) = {pt_sc+26, pt_sc+21};
+ Line(l_sc+10) = {pt_sc+28, pt_sc+23};
+ Line(l_sc+11) = {pt_sc+25, pt_sc+20};
+ Line(l_sc+12) = {pt_sc+27, pt_sc+22};
+ Line Loop(ll_sc+13) = {l_sc+8, 5 +l_sc , 6+l_sc, 7+l_sc};Plane Surface(s_sc+14) = {ll_sc+13};
+ Line Loop(ll_sc+15) = {l_sc+4, 1 +l_sc , 2+l_sc, 3+l_sc};Plane Surface(s_sc+16) = {ll_sc+15};
+ Line Loop(ll_sc+17) = {l_sc+8, 10+l_sc , -4-l_sc, -11-l_sc};Surface(s_sc+18) = {ll_sc+17};
+ Line Loop(ll_sc+19) = {l_sc+5, 9 +l_sc , -1-l_sc, -10-l_sc};Surface(s_sc+20) = {ll_sc+19};
+ Line Loop(ll_sc+21) = {l_sc+9, 2 +l_sc , -12-l_sc, -6-l_sc};Surface(s_sc+22) = {ll_sc+21};
+ Line Loop(ll_sc+23) = {l_sc+3, -11-l_sc , -7-l_sc, 12+l_sc};Surface(s_sc+24) = {ll_sc+23};
+ Surface Loop(1072) = {s_sc+14, s_sc+18, s_sc+20, s_sc+22, s_sc+16, s_sc+24};
+ Volume(1092) = {1072};
+ EndIf
+ If(flag_shape==CONE)
+ Point(pt_sc+20) = { 0, cone_ry, -cone_h/3, In_lc};
+ Point(pt_sc+21) = { 0,-cone_ry, -cone_h/3, In_lc};
+ Point(pt_sc+22) = { cone_rx, 0, -cone_h/3, In_lc};
+ Point(pt_sc+23) = {-cone_rx, 0, -cone_h/3, In_lc};
+ Point(pt_sc+24) = {0, 0, -cone_h/3, In_lc};
+ Point(pt_sc+26) = {0, 0, 2*cone_h/3, In_lc};
+ Ellipse(l_sc+1) = {pt_sc+23, pt_sc+24, pt_sc+21, pt_sc+21};
+ Ellipse(l_sc+2) = {pt_sc+21, pt_sc+24, pt_sc+22, pt_sc+22};
+ Ellipse(l_sc+3) = {pt_sc+22, pt_sc+24, pt_sc+20, pt_sc+20};
+ Ellipse(l_sc+4) = {pt_sc+20, pt_sc+24, pt_sc+24, pt_sc+23};
+ Line(l_sc+5)={pt_sc+20,pt_sc+26};Line(l_sc+6)={pt_sc+21,pt_sc+26};
+ Line(l_sc+7)={pt_sc+22,pt_sc+26};Line(l_sc+8)={pt_sc+23,pt_sc+26};
+ Line Loop(ll_sc+9) = {l_sc+1, l_sc+2, l_sc+3, l_sc+4};Plane Surface(s_sc+10) = {ll_sc+9};
+ Line Loop(ll_sc+11) = {l_sc+6, -8-l_sc, 1+l_sc};Surface(s_sc+12) = {ll_sc+11};
+ Line Loop(ll_sc+13) = {l_sc+8, -5-l_sc, 4+l_sc};Surface(s_sc+14) = {ll_sc+13};
+ Line Loop(ll_sc+15) = {l_sc+5, -7-l_sc, 3+l_sc};Surface(s_sc+16) = {ll_sc+15};
+ Line Loop(ll_sc+17) = {l_sc+7, -6-l_sc, 2+l_sc};Surface(s_sc+18) = {ll_sc+17};
+ Surface Loop(1072) = {s_sc+12, s_sc+18, s_sc+16, s_sc+14, s_sc+10};
+ Volume(1092) = {1072};
+ EndIf
+ If (flag_shape==TOR)
+ Point(pt_sc+20) = { 0, tor_r1+tor_r2x+tor_r2x , 0, In_lc};
+ Point(pt_sc+21) = { 0, tor_r1+tor_r2x-tor_r2x , 0, In_lc};
+ Point(pt_sc+22) = { 0, tor_r1+tor_r2x , tor_r2z, In_lc};
+ Point(pt_sc+23) = { 0, tor_r1+tor_r2x , -tor_r2z, In_lc};
+ Point(pt_sc+24) = { 0, tor_r1+tor_r2x, 0, In_lc};
+ Ellipse(l_sc+1) = {pt_sc+21, pt_sc+24, pt_sc+22, pt_sc+22};
+ Ellipse(l_sc+2) = {pt_sc+22, pt_sc+24, pt_sc+20, pt_sc+20};
+ Ellipse(l_sc+3) = {pt_sc+20, pt_sc+24, pt_sc+23, pt_sc+23};
+ Ellipse(l_sc+4) = {pt_sc+23, pt_sc+24, pt_sc+21, pt_sc+21};
+ Line Loop(ll_sc+5) = {l_sc+2, l_sc+3, l_sc+4, l_sc+1};Plane Surface(s_sc+6) = {ll_sc+5};
+ Extrude {{0, 0, 1}, {0, 0, 0}, deg2rad*tor_angle/2} {Surface{s_sc+6};}
+ Extrude {{0, 0, 1}, {0, 0, 0},-deg2rad*tor_angle/2} {Surface{s_sc+6};}
+ Surface Loop(1072) = {472, 460, 464, 468, 473, 494, 482, 486, 490, 495};
+ EndIf
+ inerpml_ll=newll;
+ Surface Loop(inerpml_ll) = {326, 272, 218, 362, 202, 254};
+ Volume(1093) = {inerpml_ll,1072};
+
+ Physical Point(2000) = {1}; // PrintPoint
+ Physical Volume(1000) = {402, 398, 394, 408, 442, 428, 432, 436}; // PML XYZ
+ Physical Volume(1001) = {400, 410, 430, 444}; // PML XZ
+ Physical Volume(1002) = {406, 396, 434, 440}; // PML YZ
+ Physical Volume(1003) = {412, 416, 420, 424}; // PML XY
+ Physical Volume(1004) = {404, 438}; // PML Z
+ Physical Volume(1005) = {426, 418}; // PML Y
+ Physical Volume(1006) = {414, 422}; // PML X
+ Physical Volume(1007) = {1093}; // Out
+ If (flag_shape==TOR)
+ Physical Volume(1008) = {445,446}; // Scatterer
+ Else
+ Physical Volume(1008) = {1092}; // Scatterer
+ EndIf
+
+ // Physical Volume(1007) = {392}; // Out
+ // Physical Volume(1008) = {1039,1041,1043,1045,1047,1049,1051,1053}; // Scatterer
+EndIf
+If(flag_cartpml==0)
+ Point(1) = { r_pml_in, 0, 0 , Out_lc};
+ Point(2) = {-r_pml_in, 0, 0 , Out_lc};
+ Point(3) = {0, r_pml_in, 0 , Out_lc};
+ Point(4) = {0, -r_pml_in, 0 , Out_lc};
+ Point(5) = {0, 0, r_pml_in , Out_lc};
+ Point(6) = {0, 0, -r_pml_in , Out_lc};
+ Point(7) = { r_pml_out, 0, 0, PML_lc};
+ Point(8) = {-r_pml_out, 0, 0, PML_lc};
+ Point(9) = {0, r_pml_out, 0, PML_lc};
+ Point(10) = {0, -r_pml_out, 0, PML_lc};
+ Point(11) = {0, 0, r_pml_out, PML_lc};
+ Point(12) = {0, 0, -r_pml_out, PML_lc};
+ Point(13) = { 0, 0, 0, CenterScat_lc};
+ If(flag_shape==ELL)
+ Point(20) = { 0, 0, ell_rz, In_lc}; Point(21) = { 0, 0,-ell_rz, In_lc};
+ Point(22) = { 0, ell_ry, 0, In_lc}; Point(23) = { 0,-ell_ry, 0, In_lc};
+ Point(25) = { ell_rx, 0, 0, In_lc}; Point(26) = {-ell_rx, 0, 0, In_lc};
+ Ellipse(145) = {22, 13, 26, 26} Plane{0,0,1};
+ Ellipse(146) = {26, 13, 23, 23} Plane{0,0,1};
+ Ellipse(147) = {23, 13, 25, 25} Plane{0,0,1};
+ Ellipse(148) = {25, 13, 22, 22} Plane{0,0,1};
+ Ellipse(149) = {21, 13, 26, 26} Plane{0,0,1};
+ Ellipse(150) = {26, 13, 20, 20} Plane{0,0,1};
+ Ellipse(151) = {20, 13, 25, 25} Plane{0,0,1};
+ Ellipse(152) = {25, 13, 21, 21} Plane{0,0,1};
+ Ellipse(153) = {23, 13, 20, 20} Plane{0,0,1};
+ Ellipse(154) = {20, 13, 22, 22} Plane{0,0,1};
+ Ellipse(155) = {22, 13, 21, 21} Plane{0,0,1};
+ Ellipse(156) = {21, 13, 23, 23} Plane{0,0,1};
+ Line Loop(157) = {151, -147, 153}; Surface(158) = {157};
+ Line Loop(159) = {151, 148, -154}; Surface(160) = {159};
+ Line Loop(161) = {145, 150, 154}; Surface(162) = {161};
+ Line Loop(163) = {150, -153, -146}; Surface(164) = {163};
+ Line Loop(165) = {156, 147, 152}; Surface(166) = {165};
+ Line Loop(167) = {152, -155, -148}; Surface(168) = {167};
+ Line Loop(169) = {155, 149, -145}; Surface(170) = {169};
+ Line Loop(171) = {149, 146, -156}; Surface(172) = {171};
+ Surface Loop(1072) = {160, 162, 164, 158, 166, 170, 172, 168};
+ Volume(1092) = {1072};
+ EndIf
+ If(flag_shape==PARALL)
+ Point(20) = {-par_ax/2,-par_ay/2,-par_az/2,In_lc};
+ Point(21) = { par_ax/2,-par_ay/2,-par_az/2,In_lc};
+ Point(22) = { par_ax/2, par_ay/2,-par_az/2,In_lc};
+ Point(23) = {-par_ax/2, par_ay/2,-par_az/2,In_lc};
+ Point(24) = {-par_ax/2,-par_ay/2, par_az/2,In_lc};
+ Point(25) = { par_ax/2,-par_ay/2, par_az/2,In_lc};
+ Point(26) = { par_ax/2, par_ay/2, par_az/2,In_lc};
+ Point(27) = {-par_ax/2, par_ay/2, par_az/2,In_lc};
+ Line(1) = {20, 24};
+ Line(2) = {24, 25};
+ Line(3) = {25, 21};
+ Line(4) = {21, 20};
+ Line(5) = {23, 27};
+ Line(6) = {27, 26};
+ Line(7) = {26, 22};
+ Line(8) = {22, 23};
+ Line(9) = {23, 20};
+ Line(10) = {22, 21};
+ Line(11) = {26, 25};
+ Line(12) = {27, 24};
+ Line Loop(13) = {2, -11, -6, 12};Plane Surface(14) = {13};
+ Line Loop(15) = {1, -12, -5, 9};Plane Surface(16) = {15};
+ Line Loop(17) = {3, -10, -7, 11};Plane Surface(18) = {17};
+ Line Loop(19) = {8, 5, 6, 7};Plane Surface(20) = {19};
+ Line Loop(21) = {4, -9, -8, 10};Plane Surface(22) = {21};
+ Line Loop(23) = {1, 2, 3, 4};Plane Surface(24) = {23};
+ Surface Loop(1072) = {16, 24, 14, 18, 22, 20};
+ Volume(1092) = {1072};
+ EndIf
+ If(flag_shape==CYL)
+ Point(20) = { 0, cyl_ry, -cyl_h/2, In_lc};
+ Point(21) = { 0,-cyl_ry, -cyl_h/2, In_lc};
+ Point(22) = { cyl_rx, 0, -cyl_h/2, In_lc};
+ Point(23) = {-cyl_rx, 0, -cyl_h/2, In_lc};
+ Point(24) = {0, 0, -cyl_h/2, In_lc};
+ Point(25) = { 0, cyl_ry, cyl_h/2, In_lc};
+ Point(26) = { 0,-cyl_ry, cyl_h/2, In_lc};
+ Point(27) = { cyl_rx, 0, cyl_h/2, In_lc};
+ Point(28) = {-cyl_rx, 0, cyl_h/2, In_lc};
+ Point(29) = {0, 0, cyl_h/2, In_lc};
+ Ellipse(1) = {23, 24, 21, 21};
+ Ellipse(2) = {21, 24, 22, 22};
+ Ellipse(3) = {22, 24, 20, 20};
+ Ellipse(4) = {20, 24, 24, 23};
+ Ellipse(5) = {28, 29, 26, 26};
+ Ellipse(6) = {26, 29, 27, 27};
+ Ellipse(7) = {27, 29, 25, 25};
+ Ellipse(8) = {25, 29, 28, 28};
+ Line(9) = {26, 21};
+ Line(10) = {28, 23};
+ Line(11) = {25, 20};
+ Line(12) = {27, 22};
+ Line Loop(13) = {8, 5, 6, 7};Plane Surface(14) = {13};
+ Line Loop(15) = {4, 1, 2, 3};Plane Surface(16) = {15};
+ Line Loop(17) = {8, 10, -4, -11};Surface(18) = {17};
+ Line Loop(19) = {5, 9, -1, -10};Surface(20) = {19};
+ Line Loop(21) = {9, 2, -12, -6};Surface(22) = {21};
+ Line Loop(23) = {3, -11, -7, 12};Surface(24) = {23};
+ Surface Loop(1072) = {14, 18, 20, 22, 16, 24};
+ Volume(1092) = {1072};
+ EndIf
+ If(flag_shape==CONE)
+ Point(20) = { 0, cone_ry, -cone_h/3, In_lc};
+ Point(21) = { 0,-cone_ry, -cone_h/3, In_lc};
+ Point(22) = { cone_rx, 0, -cone_h/3, In_lc};
+ Point(23) = {-cone_rx, 0, -cone_h/3, In_lc};
+ Point(24) = {0, 0, -cone_h/3, In_lc};
+ Point(26) = {0, 0, 2*cone_h/3, In_lc};
+ Ellipse(1) = {23, 24, 21, 21};
+ Ellipse(2) = {21, 24, 22, 22};
+ Ellipse(3) = {22, 24, 20, 20};
+ Ellipse(4) = {20, 24, 24, 23};
+ Line(5)={20,26};Line(6)={21,26};
+ Line(7)={22,26};Line(8)={23,26};
+ Line Loop(9) = {1, 2, 3, 4};Plane Surface(10) = {9};
+ Line Loop(11) = {6, -8, 1};Surface(12) = {11};
+ Line Loop(13) = {8, -5, 4};Surface(14) = {13};
+ Line Loop(15) = {5, -7, 3};Surface(16) = {15};
+ Line Loop(17) = {7, -6, 2};Surface(18) = {17};
+ Surface Loop(1072) = {12, 18, 16, 14, 10};
+ Volume(1092) = {1072};
+ EndIf
+ If (flag_shape==TOR)
+ Point(20) = { 0, tor_r1+tor_r2x+tor_r2x , 0, In_lc};
+ Point(21) = { 0, tor_r1+tor_r2x-tor_r2x , 0, In_lc};
+ Point(22) = { 0, tor_r1+tor_r2x , tor_r2z, In_lc};
+ Point(23) = { 0, tor_r1+tor_r2x , -tor_r2z, In_lc};
+ Point(24) = { 0, tor_r1+tor_r2x, 0, In_lc};
+ Ellipse(1) = {21, 24, 22, 22};
+ Ellipse(2) = {22, 24, 20, 20};
+ Ellipse(3) = {20, 24, 23, 23};
+ Ellipse(4) = {23, 24, 21, 21};
+ Line Loop(5) = {2, 3, 4, 1};Plane Surface(6) = {5};
+ Extrude {{0, 0, 1}, {0, 0, 0}, deg2rad*tor_angle/2} {Surface{6};}
+ Extrude {{0, 0, 1}, {0, 0, 0},-deg2rad*tor_angle/2} {Surface{6};}
+ Surface Loop(1072) = {37, 41, 45, 49, 50, 15, 19, 23, 27, 28};
+ EndIf
+ // Rotate {{0, 0, 1}, {0, 0, 0}, rot_phi*deg2rad} {Volume{1092};}
+ // Rotate {{-Sin[rot_phi*deg2rad], Cos[rot_phi*deg2rad], 0}, {0, 0, 0}, rot_theta*deg2rad} {Volume{1092};}
+
+ Circle(114) = {2 , 13, 4};
+ Circle(115) = {4 , 13, 1};
+ Circle(116) = {1 , 13, 3};
+ Circle(117) = {3 , 13, 2};
+ Circle(118) = {5 , 13, 3};
+ Circle(119) = {3 , 13, 6};
+ Circle(120) = {6 , 13, 4};
+ Circle(121) = {4 , 13, 5};
+ Circle(122) = {6 , 13, 2};
+ Circle(123) = {2 , 13, 5};
+ Circle(124) = {5 , 13, 1};
+ Circle(125) = {1 , 13, 6};
+ Circle(126) = {8 , 13, 10};
+ Circle(127) = {10, 13, 7};
+ Circle(128) = {7 , 13, 9};
+ Circle(129) = {9 , 13, 8};
+ Circle(130) = {8 , 13, 11};
+ Circle(131) = {11, 13, 7};
+ Circle(132) = {7 , 13, 12};
+ Circle(133) = {12, 13, 8};
+ Circle(134) = {11, 13, 10};
+ Circle(135) = {10, 13, 12};
+ Circle(136) = {12, 13, 9};
+ Circle(137) = {9 , 13, 11};
+
+ Line Loop(1054) = {123, 118, 117};Surface(1055) = {1054};
+ Line Loop(1056) = {118, -116, -124};Surface(1057) = {1056};
+ Line Loop(1058) = {116, 119, -125};Surface(1059) = {1058};
+ Line Loop(1061) = {119, 122, -117};Surface(1062) = {1061};
+ Line Loop(1063) = {114, -120, 122};Surface(1064) = {1063};
+ Line Loop(1065) = {114, 121, -123};Surface(1066) = {1065};
+ Line Loop(1067) = {121, 124, -115};Surface(1068) = {1067};
+ Line Loop(1069) = {125, 120, 115};Surface(1070) = {1069};
+
+ Line Loop(1074) = {128, 137, 131}; Surface(1075) = {1074};
+ Line Loop(1076) = {129, 130, -137}; Surface(1077) = {1076};
+ Line Loop(1078) = {136, 129, -133}; Surface(1079) = {1078};
+ Line Loop(1080) = {128, -136, -132}; Surface(1081) = {1080};
+ Line Loop(1082) = {134, 127, -131}; Surface(1083) = {1082};
+ Line Loop(1084) = {127, 132, -135}; Surface(1085) = {1084};
+ Line Loop(1086) = {135, 133, 126}; Surface(1087) = {1086};
+ Line Loop(1088) = {126, -134, -130}; Surface(1089) = {1088};
+ // pml in
+ Surface Loop(1100) = {1057, 1055, 1066, 1064, 1070, 1059, 1062, 1068};
+ // pml out
+ Surface Loop(1101) = {1089, 1087, 1085, 1083, 1075, 1081, 1079, 1077};
+
+
+ Surface Loop(1071) = {1062, 1059, 1057, 1068, 1066, 1064, 1070, 1055};
+ Volume(1073) = {1071, 1072};
+ Surface Loop(1090) = {1075, 1081, 1079, 1077, 1089, 1087, 1085, 1083};
+ Volume(1091) = {1071, 1090};
+ If (flag_shape==TOR)
+ Physical Volume(1) = {1,2}; // Scatterer
+ Else
+ Physical Volume(1) = {1092}; // Scatterer
+ EndIf
+ Physical Volume(2) = {1073}; // Out - air
+ Physical Volume(3) = {1091}; // pml
+ Physical Point(2000) = {1}; // PrintPoint
+EndIf
+
+// Mesh.Algorithm = 1; // // 1=MeshAdapt, 5=Delaunay, 6=Frontal
+// Mesh.Algorithm3D = 1;//4; // // 1=Delaunay, 4=Frontal
+Mesh.Optimize = 1;
+Solver.AutoMesh=2;
+Geometry.Points = 1;
+// Mesh.SurfaceEdges = 0;
+Mesh.VolumeEdges = 0;
+
+Printf(Sprintf("nm = %.9e;",nm)) > "scattering_tmp.py";
+Printf(Sprintf("lambda0 = %.9e;",lambda0)) >> "scattering_tmp.py";
+Printf(Sprintf("rbb = %.9e;",rbb)) >> "scattering_tmp.py";
+Printf(Sprintf("epsr_In_re = %.9e;",epsr_In_re)) >> "scattering_tmp.py";
+Printf(Sprintf("epsr_In_im = %.9e;",epsr_In_im)) >> "scattering_tmp.py";
+Printf(Sprintf("epsr_Out_re = %.9e;",epsr_Out_re)) >> "scattering_tmp.py";
+Printf(Sprintf("epsr_Out_im = %.9e;",epsr_Out_im)) >> "scattering_tmp.py";
+Printf(Sprintf("sph_scan = %.9e;",sph_scan)) >> "scattering_tmp.py";
+Printf(Sprintf("r_sph_min = %.9e;",r_sph_min)) >> "scattering_tmp.py";
+Printf(Sprintf("r_sph_max = %.9e;",r_sph_max)) >> "scattering_tmp.py";
+Printf(Sprintf("nb_cuts = %g;" ,nb_cuts)) >> "scattering_tmp.py";
+Printf(Sprintf("npts_theta = %g;" ,npts_theta)) >> "scattering_tmp.py";
+Printf(Sprintf("npts_phi = %g;" ,npts_phi)) >> "scattering_tmp.py";
+Printf(Sprintf("flag_study = %g;" ,flag_study)) >> "scattering_tmp.py";
+Printf(Sprintf("n_max = %g;" ,n_max)) >> "scattering_tmp.py";
+Printf(Sprintf("p_max = %g;" ,p_max)) >> "scattering_tmp.py";
+Printf(Sprintf("siwt = %g;" ,siwt)) >> "scattering_tmp.py";
+
diff --git a/ElectromagneticScattering/scattering.pro b/ElectromagneticScattering/scattering.pro
new file mode 100644
index 0000000000000000000000000000000000000000..98159da821037281a6b4f75d9955f13db5d9f2ee
--- /dev/null
+++ b/ElectromagneticScattering/scattering.pro
@@ -0,0 +1,659 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering.pro //
+///////////////////////////////
+
+Include "scattering_data.geo";
+myDir = "run_results/";
+
+Group {
+ If(flag_cartpml==1)
+ PMLxyz = Region[1000];
+ PMLxz = Region[1001];
+ PMLyz = Region[1002];
+ PMLxy = Region[1003];
+ PMLz = Region[1004];
+ PMLy = Region[1005];
+ PMLx = Region[1006];
+ Scat_In = Region[1008];
+ Scat_Out = Region[1007];
+ // Domains
+ Domain = Region[{Scat_In,Scat_Out}];
+ PMLs = Region[{PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
+ All_domains = Region[{Scat_In,Scat_Out,PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
+ EndIf
+ If(flag_cartpml==0)
+ Scat_In = Region[1];
+ Scat_Out = Region[2];
+ Domain = Region[{Scat_In,Scat_Out}];
+ PMLs = Region[3];
+ All_domains = Region[{Scat_In,Scat_Out,PMLs}];
+ EndIf
+ PrintPoint = Region[2000];
+}
+
+Function{
+ // test={1.4,2.6,3.7};
+ // For i In {0:#test()-1}
+ // Printf("%g",test(i));
+ // EndFor
+ I[] = Complex[0,1];
+ avoid_sing = 1.e-14;
+ mu0 = 4*Pi*100.0*nm;
+ epsilon0 = 8.854187817e-3*nm;
+ cel = 1.0/Sqrt[epsilon0 * mu0];
+ Freq = cel/lambda0;
+ omega0 = 2.0*Pi*Freq;
+ k_Out = 2.0*Pi*Sqrt[epsr_Out_re]/lambda0;
+ r3D_sph[] = Sqrt[X[]^2+Y[]^2+Z[]^2];
+ r2D_sph[] = Sqrt[X[]^2+Y[]^2];
+ cos_theta[] = Z[]/(r3D_sph[]+avoid_sing);
+ theta[] = Acos[cos_theta[]];
+ phi[] = Atan2[-Y[],-X[]]+Pi;
+
+ For kr In {0:nb_cuts-1}
+ r_sph~{kr}=r_sph_min+kr*(r_sph_max-r_sph_min)/(nb_cuts-1);
+ EndFor
+
+ If (flag_study==0)
+ Ae = 1.0;
+ Ah = Ae*Sqrt[epsilon0*epsr_In_re/mu0];
+ alpha0 = k_Out*Sin[theta0]*Cos[phi0];
+ beta0 = k_Out*Sin[theta0]*Sin[phi0];
+ gamma0 = k_Out*Cos[theta0];
+ Ex0 = Ae * Cos[psi0]*Cos[theta0]*Cos[phi0] - Ae* Sin[psi0]*Sin[phi0];
+ Ey0 = Ae * Cos[psi0]*Cos[theta0]*Sin[phi0] + Ae* Sin[psi0]*Cos[phi0];
+ Ez0 = -Ae * Cos[psi0]*Sin[theta0];
+ Hx0 = -1/(omega0*mu0)*(beta0 * Ez0 - gamma0 * Ey0);
+ Hy0 = -1/(omega0*mu0)*(gamma0 * Ex0 - alpha0 * Ez0);
+ Hz0 = -1/(omega0*mu0)*(alpha0 * Ey0 - beta0 * Ex0);
+ Prop[] = Ae * Complex[ Cos[alpha0*X[]+beta0*Y[]+gamma0*Z[]] , siwt*Sin[alpha0*X[]+beta0*Y[]+gamma0*Z[]] ];
+ Einc[] = Vector[Ex0*Prop[],Ey0*Prop[],Ez0*Prop[]];
+ Hinc[] = Vector[Hx0*Prop[],Hy0*Prop[],Hz0*Prop[]];
+ Pinc = 0.5*Ae*Ae*Sqrt[epsilon0*epsr_In_re/mu0] * Cos[theta0];
+ EndIf
+ If (flag_study==2)
+ // Green Dyadic
+ normrr_p[] = Sqrt[(X[]-x_p)^2+(Y[]-y_p)^2+(Z[]-z_p)^2];
+ Gsca[] = Complex[Cos[-1.*siwt*k_Out*normrr_p[]],
+ Sin[-1.*siwt*k_Out*normrr_p[]]]
+ /(4.*Pi*normrr_p[]);
+ GDfact_diag[] = 1.+(I[]*k_Out*normrr_p[]-1.)/(k_Out*normrr_p[])^2;
+ GDfact_glob[] = (3.-3.*I[]*k_Out*normrr_p[]-(k_Out*normrr_p[])^2)/
+ (k_Out*normrr_p[]^2)^2;
+ Dyad_Green2[]=
+ Gsca[]*(
+ GDfact_diag[]*TensorDiag[1,1,1]
+ + GDfact_glob[]*
+ Tensor[(X[]-x_p)*(X[]-x_p),(X[]-x_p)*(Y[]-y_p),(X[]-x_p)*(Z[]-z_p),
+ (Y[]-y_p)*(X[]-x_p),(Y[]-y_p)*(Y[]-y_p),(Y[]-y_p)*(Z[]-z_p),
+ (Z[]-z_p)*(X[]-x_p),(Z[]-z_p)*(Y[]-y_p),(Z[]-z_p)*(Z[]-z_p)]);
+ // Getdp Func
+ Dyad_Green[] = DyadGreenHom[x_p,y_p,z_p,XYZ[],k_Out,siwt];
+ Einc_Green_0[] = Dyad_Green2[]*Vector[1,0,0];
+ Einc_Green_1[] = Dyad_Green2[]*Vector[0,1,0];
+ Einc_Green_2[] = Dyad_Green2[]*Vector[0,0,1];
+ CurlDyad_Green[] = CurlDyadGreenHom[x_p,y_p,z_p,XYZ[],k_Out,siwt];
+ Hinc_Green_0[] = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[1,0,0];
+ Hinc_Green_1[] = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[0,1,0];
+ Hinc_Green_2[] = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[0,0,1];
+ // test[] = DyadGreenHom[x_p,y_p,z_p,x_p+1e-1*nm,y_p+1e-1*nm,z_p+1e-1*nm,lambda0,Sqrt[epsr_Out_re]]*Vector[1,0,0];
+ EndIf
+
+ a_pml = 1.;
+ b_pml = -siwt;
+
+ eigtarget = (2.*Pi*cel/lambda0)^2;
+ eigfilter = 1e-4*Sqrt[eigtarget];
+ EigFilter[] = (Norm[$EigenvalueReal] > eigfilter);
+
+ If(flag_cartpml==1)
+ sx[Scat_In] = 1.;
+ sy[Scat_In] = 1.;
+ sz[Scat_In] = 1.;
+ sx[Scat_Out] = 1.;
+ sy[Scat_Out] = 1.;
+ sz[Scat_Out] = 1.;
+ sx[PMLxyz] = Complex[a_pml,b_pml];
+ sy[PMLxyz] = Complex[a_pml,b_pml];
+ sz[PMLxyz] = Complex[a_pml,b_pml];
+ sx[PMLxz] = Complex[a_pml,b_pml];
+ sy[PMLxz] = 1.0;
+ sz[PMLxz] = Complex[a_pml,b_pml];
+ sx[PMLyz] = 1.0;
+ sy[PMLyz] = Complex[a_pml,b_pml];
+ sz[PMLyz] = Complex[a_pml,b_pml];
+ sx[PMLxy] = Complex[a_pml,b_pml];
+ sy[PMLxy] = Complex[a_pml,b_pml];
+ sz[PMLxy] = 1.0;
+ sx[PMLx] = Complex[a_pml,b_pml];
+ sy[PMLx] = 1.0;
+ sz[PMLx] = 1.0;
+ sx[PMLy] = 1.0;
+ sy[PMLy] = Complex[a_pml,b_pml];
+ sz[PMLy] = 1.0;
+ sx[PMLz] = 1.0;
+ sy[PMLz] = 1.0;
+ sz[PMLz] = Complex[a_pml,b_pml];
+ Lxx[] = sy[]*sz[]/sx[];
+ Lyy[] = sz[]*sx[]/sy[];
+ Lzz[] = sx[]*sy[]/sz[];
+
+ epsilonr_In[] = Complex[epsr_In_re , epsr_In_im];
+ epsilonr_Out[] = Complex[epsr_Out_re , epsr_Out_im];
+
+ epsilonr[Scat_In] = epsilonr_In[] * TensorDiag[1.,1.,1.];
+ epsilonr[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr[PMLs] = epsr_Out_re * TensorDiag[Lxx[],Lyy[],Lzz[]];
+
+ epsilonr1[Scat_In] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr1[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr1[PMLs] = epsr_Out_re * TensorDiag[Lxx[],Lyy[],Lzz[]];
+
+ mur[Scat_In] = TensorDiag[1.,1.,1.];
+ mur[Scat_Out] = TensorDiag[1.,1.,1.];
+ mur[PMLs] = TensorDiag[Lxx[],Lyy[],Lzz[]];
+ EndIf
+
+ If(flag_cartpml==0)
+ sr[] = Complex[a_pml,b_pml];
+ sphi[] = Complex[1,0];
+ stheta[] = Complex[1,0];
+ r_tild[] = r_pml_in + (r3D_sph[] - r_pml_in) * sr[];
+ theta_tild[] = theta[];
+ pml_tens_temp1[] = TensorDiag[(r_tild[]/r3D_sph[])^2 * sphi[]*stheta[]/sr[],
+ (r_tild[]/r3D_sph[]) * sr[]*stheta[]/sphi[],
+ (r_tild[]/r3D_sph[]) * sphi[]*sr[]/stheta[]];
+ pml_tens_temp2[] = Rotate[pml_tens_temp1[],0,-theta[]-Pi/2,0];
+ pml_tens[] = Rotate[pml_tens_temp2[],0,0,-phi[]];
+
+ epsilonr_In[] = Complex[epsr_In_re , epsr_In_im];
+ epsilonr_Out[] = Complex[epsr_Out_re , epsr_Out_im];
+
+ epsilonr[Scat_In] = epsilonr_In[] * TensorDiag[1.,1.,1.];
+ epsilonr[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr[PMLs] = pml_tens[];
+
+ epsilonr1[Scat_In] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr1[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr1[PMLs] = pml_tens[];
+
+ mur[Scat_In] = TensorDiag[1.,1.,1.];
+ mur[Scat_Out] = TensorDiag[1.,1.,1.];
+ mur[PMLs] = pml_tens[];
+ EndIf
+
+ If(flag_study==RES_PW)
+ source[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc[];
+ EndIf
+ If(flag_study==RES_TMAT)
+ For pe In {1:p_max}
+ ne = Floor[Sqrt[pe]];
+ me = ne*(ne+1) - Floor[pe];
+ Mnm_source~{pe}[] = Mnm[1,ne,me,XYZ[],k_Out];
+ Nnm_source~{pe}[] = Nnm[1,ne,me,XYZ[],k_Out];
+ source_M~{pe}[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Mnm_source~{pe}[];
+ source_N~{pe}[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm_source~{pe}[];
+ EndFor
+ EndIf
+
+ If (flag_study==RES_GREEN)
+ sourceG_0[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_0[];
+ sourceG_1[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_1[];
+ sourceG_2[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_2[];
+ EndIf
+}
+
+Constraint {
+ // {Name Dirichlet; Type Assign;
+ // Case {
+ // { Region SurfDirichlet; Value 0.; }
+ // }
+ // }
+}
+
+Jacobian {
+ { Name JVol ;
+ Case {
+ { Region All ; Jacobian Vol ; }
+ }
+ }
+ { Name JSur ;
+ Case {
+ { Region All ; Jacobian Sur ; }
+ }
+ }
+ { Name JLin ;
+ Case {
+ { Region All ; Jacobian Lin ; }
+ }
+ }
+}
+
+Integration {
+ { Name Int_1 ;
+ Case {
+ { Type Gauss ;
+ Case {
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line ; NumberOfPoints 4 ; }
+ { GeoElement Triangle ; NumberOfPoints 6 ; }
+ { GeoElement Tetrahedron ; NumberOfPoints 15 ; }
+ }
+ }
+ }
+ }
+}
+
+FunctionSpace {
+ { Name Hcurl; Type Form1;
+ BasisFunction {
+ { Name sn; NameOfCoef un; Function BF_Edge;
+ Support Region[All_domains]; Entity EdgesOf[All]; }
+ { Name sn2; NameOfCoef un2; Function BF_Edge_2E;
+ Support Region[All_domains]; Entity EdgesOf[All]; }
+ If (is_FEM_o2==1)
+ { Name sn3; NameOfCoef un3; Function BF_Edge_3F_b;
+ Support Region[All_domains]; Entity FacetsOf[All]; }
+ { Name sn4; NameOfCoef un4; Function BF_Edge_3F_c;
+ Support Region[All_domains]; Entity FacetsOf[All]; }
+ { Name sn5; NameOfCoef un5; Function BF_Edge_4E;
+ Support Region[All_domains]; Entity EdgesOf[All]; }
+ EndIf
+ }
+ Constraint {
+ // { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ // { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ }
+ }
+}
+
+Formulation {
+ If (flag_study==RES_QNM)
+ {Name QNM_helmholtz_vector; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [ -1/mur[] * Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { DtDtDof[-1. * 1/cel^2 * epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ EndIf
+
+ If (flag_study==RES_PW)
+ {Name PW_helmholtz_vector; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [source[] , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ EndIf
+ If (flag_study==RES_TMAT)
+ For pe In {1:p_max}
+ {Name VPWM_helmholtz_vector~{pe}; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [source_M~{pe}[] , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ {Name VPWN_helmholtz_vector~{pe}; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [source_N~{pe}[] , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ EndFor
+ EndIf
+ If (flag_study==RES_GREEN)
+ For ncomp In {0:2}
+ {Name GreenAll_helmholtz_vector~{ncomp}; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [sourceG~{ncomp}[] , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ EndFor
+ EndIf
+}
+
+Resolution {
+ If (flag_study==RES_PW)
+ { Name res_PW_helmholtz_vector;
+ System {
+ { Name P; NameOfFormulation PW_helmholtz_vector; Type ComplexValue; Frequency Freq; }
+ }
+ Operation {
+ CreateDir[Str[myDir]];
+ Evaluate[Python[]{"scattering_init.py"}];
+ Generate[P];
+ Solve[P];
+ PostOperation[PW_postop];
+ Evaluate[Python[]{"scattering_post.py"}];
+ }
+ }
+ EndIf
+ If (flag_study==RES_TMAT)
+ { Name res_VPWall_helmholtz_vector;
+ System {
+ For pe In {1:p_max}
+ { Name M~{pe}; NameOfFormulation VPWM_helmholtz_vector~{pe}; Type ComplexValue; Frequency Freq; }
+ { Name N~{pe}; NameOfFormulation VPWN_helmholtz_vector~{pe}; Type ComplexValue; Frequency Freq; }
+ EndFor
+ }
+ Operation {
+ CreateDir[Str[myDir]];
+ Evaluate[Python[]{"scattering_init.py"}];
+ For pe In {1:p_max}
+ Generate[M~{pe}];
+ Solve[M~{pe}];
+ PostOperation[VPWM_postop~{pe}];
+ Generate[N~{pe}];
+ Solve[N~{pe}];
+ PostOperation[VPWN_postop~{pe}];
+ EndFor
+ Evaluate[Python[]{"scattering_post.py"}];
+ }
+ }
+ EndIf
+ If (flag_study==RES_GREEN)
+ { Name res_GreenAll_helmholtz_vector;
+ System {
+ For ncomp In {0:2}
+ { Name M~{ncomp}; NameOfFormulation GreenAll_helmholtz_vector~{ncomp}; Type ComplexValue; Frequency Freq; }
+ EndFor
+ }
+ Operation {
+ CreateDir[Str[myDir]];
+ Evaluate[Python[]{"scattering_init.py"}];
+ For ncomp In {0:2}
+ Generate[M~{ncomp}];
+ Solve[M~{ncomp}];
+ PostOperation[GreenAll_postop~{ncomp}];
+ EndFor
+ Evaluate[Python[]{"scattering_post.py"}];
+ }
+ }
+ EndIf
+ If (flag_study==RES_QNM)
+ { Name res_QNM_helmholtz_vector;
+ System{
+ { Name E; NameOfFormulation QNM_helmholtz_vector; Type ComplexValue; }
+ }
+ Operation{
+ CreateDir[Str[myDir]];
+ Evaluate[Python[]{"scattering_init.py"}];
+ GenerateSeparate[E];
+ EigenSolve[E,neig,eigtarget,0,EigFilter[]];
+ PostOperation[QNM_postop];
+ Evaluate[Python[]{"scattering_post.py"}];
+ }
+ }
+ EndIf
+
+}
+
+PostProcessing {
+ If (flag_study==RES_QNM)
+ { Name QNM_postpro; NameOfFormulation QNM_helmholtz_vector; NameOfSystem E;
+ Quantity {
+ { Name u ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JVol; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JVol; } } }
+ }
+ }
+ EndIf
+
+ If (flag_study==RES_PW)
+ { Name PW_postpro ; NameOfFormulation PW_helmholtz_vector;NameOfSystem P;
+ Quantity {
+ { Name E_tot ; Value { Local { [{u}+Einc[]]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat_sph ; Value { Local { [Vector[
+ (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+ ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name H_scat ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}] ; In All_domains; Jacobian JVol; } } }
+ { Name H_tot ; Value { Local { [ siwt*I[]/(mur[]*mu0*omega0)*{Curl u} +Hinc[]]; In All_domains; Jacobian JVol; } } }
+ { Name normalized_losses ; Value { Integral { [ 0.5*omega0*epsilon0*Fabs[Im[epsilonr_In[]]]*(SquNorm[{u}+Einc[]]) ] ; In Scat_In ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Poy_dif ; Value { Local { [ 0.5*Re[Cross[{u} , Conj[ siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]]] ]; In All_domains; Jacobian JVol; } } }
+ { Name Poy_tot ; Value { Local { [ 0.5*Re[Cross[{u}+Einc[] , Conj[Hinc[]+siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]]] ]; In All_domains; Jacobian JVol; } } }
+ }
+ }
+ EndIf
+ If (flag_study==RES_TMAT)
+ For pe In {1:p_max}
+ { Name VPWM_postpro~{pe}; NameOfFormulation VPWM_helmholtz_vector~{pe};NameOfSystem M~{pe};
+ Quantity {
+ { Name E_scat ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat_sph ; Value { Local { [Vector[
+ (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+ ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name H_scat ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
+ { Name Mnm_source~{pe} ; Value { Local { [ Mnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
+ }
+ }
+ { Name VPWN_postpro~{pe}; NameOfFormulation VPWN_helmholtz_vector~{pe};NameOfSystem N~{pe};
+ Quantity {
+ { Name E_scat ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat_sph ; Value { Local { [Vector[
+ (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+ ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name H_scat ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
+ { Name Nnm_source~{pe} ; Value { Local { [ Nnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
+ }
+ }
+ EndFor
+ EndIf
+ If (flag_study==RES_GREEN)
+ For ncomp In {0:2}
+ { Name GreenAll_postpro~{ncomp}; NameOfFormulation GreenAll_helmholtz_vector~{ncomp};NameOfSystem M~{ncomp};
+ Quantity {
+ { Name E_tot~{ncomp} ; Value { Local {[{u}+Einc_Green~{ncomp}[]]; In All_domains; Jacobian JVol; } } }
+ { Name H_tot~{ncomp} ; Value { Local {[siwt*I[]/(mur[]*mu0*omega0)*{Curl u}+Hinc_Green~{ncomp}[]]; In All_domains; Jacobian JVol; } } }
+ { Name Einc_Green~{ncomp}; Value { Local {[ Einc_Green~{ncomp}[] ]; In All_domains; Jacobian JVol; } } }
+ { Name Hinc_Green~{ncomp}; Value { Local {[ Hinc_Green~{ncomp}[] ]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat~{ncomp} ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name E_tot_im~{ncomp} ; Value { Local {[Im[{u}+Einc_Green~{ncomp}[]]]; In All_domains; Jacobian JVol; } } }
+ { Name Einc_Green_im~{ncomp}; Value { Local {[Im[ Einc_Green~{ncomp}[] ]]; In All_domains; Jacobian JVol; } } }
+ // { Name Einc_Green_im~{ncomp}; Value { Local {[test[]]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat_im~{ncomp} ; Value { Local {[Im[{u} ]]; In All_domains; Jacobian JVol; } } }
+ { Name E_tot_sph~{ncomp}; Value { Local { [Vector[
+ (X[]*CompX[{u}+Einc_Green~{ncomp}[]]
+ + Y[]*CompY[{u}+Einc_Green~{ncomp}[]]
+ + Z[]*CompZ[{u}+Einc_Green~{ncomp}[]])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}+Einc_Green~{ncomp}[]]
+ + Y[]*Z[]*CompY[{u}+Einc_Green~{ncomp}[]]
+ - (X[]^2+Y[]^2)*CompZ[{u}+Einc_Green~{ncomp}[]])
+ /(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}+Einc_Green~{ncomp}[]]
+ + X[]*CompY[{u}+Einc_Green~{ncomp}[]])
+ /r2D_sph[] ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name E_scat_sph~{ncomp}; Value { Local { [Vector[
+ (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}]
+ - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+ ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name H_scat~{ncomp} ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
+ }
+ }
+ EndFor
+ EndIf
+}
+
+
+
+PostOperation {
+ If (flag_study==RES_QNM)
+ { Name QNM_postop; NameOfPostProcessing QNM_postpro ;
+ Operation {
+ Print [EigenValuesReal, OnElementsOf PrintPoint , Format TimeTable, File StrCat[myDir,"EigenValuesReal.txt"]];
+ Print [EigenValuesImag, OnElementsOf PrintPoint , Format TimeTable, File StrCat[myDir,"EigenValuesImag.txt"]];
+ Print [ u , OnElementsOf All_domains, File StrCat[myDir,"eigenVectors.pos"], EigenvalueLegend];
+ }
+ }
+ EndIf
+ If (flag_study==RES_PW)
+ {Name PW_postop; NameOfPostProcessing PW_postpro ;
+ Operation {
+ Print [ E_tot , OnElementsOf Domain , File StrCat[myDir,"E_tot.pos"]];
+ Print [ E_scat , OnElementsOf All_domains, File StrCat[myDir,"E_scat.pos"]];
+ // Print [ E_scat_sph , OnElementsOf All_domains, File StrCat[myDir,"E_scat_sph.pos"]];
+ Print [ H_scat , OnElementsOf Domain , File StrCat[myDir,"H_scat.pos"]];
+ Print [ H_tot , OnElementsOf All_domains, File StrCat[myDir,"H_tot.pos"]];
+ Print [ Poy_dif , OnElementsOf All_domains, File StrCat[myDir,"Poy_dif.pos"]];
+ Print [ Poy_tot , OnElementsOf Domain , File StrCat[myDir,"Poy_tot.pos"]];
+ For kr In {0:nb_cuts-1}
+ Print [ E_scat_sph , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ File StrCat[myDir,StrCat["E_scat_onsphere_sph_PW",Sprintf["_r%g.dat",kr]]], Format Table];
+ EndFor
+ }
+ }
+ EndIf
+ If (flag_study==RES_TMAT)
+ For pe In {1:p_max}
+ {Name VPWM_postop~{pe}; NameOfPostProcessing VPWM_postpro~{pe} ;
+ Operation {
+ For kr In {0:nb_cuts-1}
+ Print [ E_scat_sph , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_sph_M",Sprintf["%g",pe]],Sprintf["_r%g.dat",kr]]], Format Table];
+ EndFor
+ Print [ E_scat , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} },
+ File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]],Sprintf["_r%g.pos",kr]]],
+ Name StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]]];
+ }
+ }
+ {Name VPWN_postop~{pe}; NameOfPostProcessing VPWN_postpro~{pe} ;
+ Operation {
+ For kr In {0:nb_cuts-1}
+ Print [ E_scat_sph , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_sph_N",Sprintf["%g",pe]],Sprintf["_r%g.dat",kr]]], Format Table];
+ EndFor
+ Print [ E_scat , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} },
+ File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_cart_N",Sprintf["%g",pe]],Sprintf["_r%g.pos",kr]]],
+ Name StrCat["E_scat_onsphere_cart_N",Sprintf["%g",pe]]];
+ }
+ }
+ EndFor
+ EndIf
+ If (flag_study==RES_GREEN)
+ For ncomp In {0:2}
+ {Name GreenAll_postop~{ncomp}; NameOfPostProcessing GreenAll_postpro~{ncomp} ;
+ Operation {
+ Print [ E_tot_sph~{ncomp} , OnGrid
+ {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+ {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ File StrCat[myDir,StrCat["E_tot_onsphere_sph_G",Sprintf["%g.dat",ncomp]]] , Format Table];
+ Print [ E_tot~{ncomp} , OnGrid
+ {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+ {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} },
+ File StrCat[myDir,StrCat["E_tot_onsphere_cart_G",Sprintf["%g.pos",ncomp]]]];
+ Print [ E_scat_im~{ncomp} , OnPoint {x_p,y_p,z_p},
+ File StrCat[myDir,StrCat["E_scat_im_onpt_G",Sprintf["%g.dat",ncomp]]] , Format Table];
+ // Print [ Einc_Green_im~{ncomp} , OnPoint {x_p+0.1*nm,y_p+0.1*nm,z_p+0.1*nm},
+ // File StrCat[myDir,StrCat["E_inc_im_onpt_G",Sprintf["%g.dat",ncomp]]] , Format Table];
+ // If(flag_FF==1)
+ // Print [ E_tot~{ncomp}, OnElementsOf Scat_Out , File StrCat[myDir,StrCat["E_tot_",Sprintf["%g.pos",ncomp]]]];
+ // Print [ H_tot~{ncomp}, OnElementsOf Scat_Out , File StrCat[myDir,StrCat["H_tot_",Sprintf["%g.pos",ncomp]]]];
+ // Print [ E_tot~{ncomp} , OnGrid
+ // {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+ // {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ // {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ // File StrCat[myDir,StrCat["E_tot_onsphere_G",Sprintf["%g.pos",ncomp]]]];
+ // Print [ H_tot~{ncomp} , OnGrid
+ // {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+ // {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ // {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ // File StrCat[myDir,StrCat["H_tot_onsphere_G",Sprintf["%g.pos",ncomp]]]];
+ //
+ // Echo[
+ // Str["Plugin(NearToFarField).Wavenumber = 2*Pi/lambda0;",
+ // "Plugin(NearToFarField).PhiStart = 0;",
+ // "Plugin(NearToFarField).PhiEnd = 2*Pi;",
+ // "Plugin(NearToFarField).NumPointsPhi = npts_phi;",
+ // "Plugin(NearToFarField).ThetaStart = 0;",
+ // "Plugin(NearToFarField).ThetaEnd = Pi;",
+ // "Plugin(NearToFarField).NumPointsTheta = npts_theta;",
+ // "Plugin(NearToFarField).EView = PostProcessing.NbViews-2;",
+ // "Plugin(NearToFarField).HView = PostProcessing.NbViews-1;",
+ // "Plugin(NearToFarField).Normalize = 0;",
+ // "Plugin(NearToFarField).dB = 0;",
+ // "Plugin(NearToFarField).NegativeTime = 0;",
+ // "Plugin(NearToFarField).RFar = r_pml_in;",
+ // "Plugin(NearToFarField).Run;"], File StrCat[myDir,"tmp1.geo"]] ;
+ // EndIf
+ }
+ }
+ EndFor
+ EndIf
+}
+
+If (flag_study==RES_QNM)
+ DefineConstant[
+ R_ = {"res_QNM_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {"-solve -pos -petsc_prealloc 500 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps", Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
+
+If (flag_study==RES_TMAT)
+ DefineConstant[
+ R_ = {"res_VPWall_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {"-solve -pos -petsc_prealloc 500 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps", Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
+If (flag_study==RES_PW)
+ DefineConstant[
+ R_ = {"res_PW_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {"-solve -pos -petsc_prealloc 500 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps", Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
+If (flag_study==RES_GREEN)
+ DefineConstant[
+ R_ = {"res_GreenAll_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {"-solve -pos -petsc_prealloc 500 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps", Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
\ No newline at end of file
diff --git a/ElectromagneticScattering/scattering_data.geo b/ElectromagneticScattering/scattering_data.geo
new file mode 100644
index 0000000000000000000000000000000000000000..91e3c08c3ced34e22ad336208c257be16360a813
--- /dev/null
+++ b/ElectromagneticScattering/scattering_data.geo
@@ -0,0 +1,215 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering_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 = "2Study 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,GmshOption "Reset", Autocheck 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 = {300 , Name StrCat[pp0 , "1cylinder X-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ cyl_ry = {300 , Name StrCat[pp0 , "2cylinder Y-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ cyl_h = {300 , 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}];
+
+DefineConstant[
+ flag_study = { RES_TMAT , Name StrCat[pp1, "0Study type"],
+ Choices {RES_PW="Plane Wave Response",
+ RES_TMAT="T-Matrix",
+ RES_GREEN="Green's Tensor",
+ RES_QNM="Quasi-Normal Modes"},Closed 0,GmshOption "Reset", Autocheck 0}];
+
+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}
+];
+
+If (flag_study!=RES_QNM)
+ DefineConstant[
+ lambda0 = {587.6 , Name StrCat[pp2 , "4wavelength [nm]"] , Highlight Str[colorppOK] , Closed 0,GmshOption "Reset", Autocheck 0}];
+Else
+ DefineConstant[
+ lambda0 = {1000 , Name StrCat[pp2 , "4wavelength target [nm]"] , Highlight Str[colorppOK] , Closed 0,GmshOption "Reset", Autocheck 0},
+ neig = {30 , Name StrCat[pp2 , "5number of eigenvalues []"] , Highlight Str[colorppOK] , Closed 0,GmshOption "Reset", Autocheck 0}
+ ];
+EndIf
+lambda_bg = lambda0/Sqrt[epsr_Out_re];
+
+If (flag_study==RES_PW)
+ DefineConstant[
+ theta0 = {0 , Name StrCat[pp2 , "5plane wave theta [deg]"], Highlight Str[colorppOK] , Closed 0, Min 0, Max 180},
+ phi0 = {0 , Name StrCat[pp2 , "6plane wave phi [deg]"] , Highlight Str[colorppOK] , Closed 0, Min 0, Max 360},
+ psi0 = {0 , Name StrCat[pp2 , "7plane wave psi [deg]"] , Highlight Str[colorppOK] , Closed 0, Min 0, Max 180}];
+ theta0 = theta0 *deg2rad;
+ phi0 = phi0 * deg2rad;
+ psi0 = psi0 * deg2rad;
+EndIf
+If (flag_study==RES_GREEN)
+ DefineConstant[
+ x_p = {0 , Name StrCat[pp2 , "5Green X-point [nm]"] , Highlight Str[colorppOK] , Closed 0},
+ y_p = {0 , Name StrCat[pp2 , "6Green Y-point [nm]"] , Highlight Str[colorppOK] , Closed 0},
+ z_p = {-rbb*1.1 , Name StrCat[pp2 , "7Green Z-point [nm]"] , Highlight Str[colorppOK] , Closed 0}];
+ x_p=x_p*nm;
+ y_p=y_p*nm;
+ z_p=z_p*nm;
+EndIf
+
+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}
+];
+DefineConstant[
+ flag_cartpml = {0 , Name StrCat[pp3 , "0PML type"] , Choices {0="spherical",1="cartesian"}, Closed 0},
+ pml_size = {lambda_bg/1.5, Name StrCat[pp3 , "1PML thickness [nm]"] , Highlight Str[colorppWA] , Closed 0, Min (lambda_bg/10), Max (3*lambda_bg)},
+ paramaille = {4. , Name StrCat[pp3 , "2mesh size"], Highlight Str[colorppWA] , Closed 0},
+ refine_scat = {2. , 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[
+ space2pml = {lambda_bg/10.*4 , Name StrCat[pp4 , "0space around scatterer [nm]"] , Highlight Str[colorppNO] , Closed 1,Min (lambda_bg/10.), Max (3*lambda_bg)},
+ dist_rcut = {lambda_bg/10. , Name StrCat[pp4 , "1min dist from object & PML"] , Highlight Str[colorppNO] , Closed 1},
+ nb_cuts = {3 , Name StrCat[pp4 , "2number of cuts [integer]"] , Highlight Str[colorppNO] , Closed 1},
+ npts_theta = {50 , Name StrCat[pp4 , "3polar sampling [integer]"], Highlight Str[colorppNO] , Closed 0,Min 50, Max 300},
+ npts_phi = {100 , Name StrCat[pp4 , "4azimuthal sampling [integer]"], Highlight Str[colorppNO] , Closed 0,Min 100, Max 600}
+];
+DefineConstant[
+ flag_plotcuts = {0, Choices{0,1}, Name StrCat[pp5, "Plot radial cuts?"]},
+ flag_FF = {1, Choices{0,1}, Name StrCat[pp5, "Plot far field?"]}
+];
+
+
+
+// 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;
+dist_rcut = dist_rcut*nm;
+space2pml = space2pml*nm;
+rbb = rbb*nm;
+
+r_pml_in = space2pml+rbb;
+r_pml_out = space2pml+rbb+pml_size;
+r_sph_min = rbb+dist_rcut;
+r_sph_max = space2pml+rbb-dist_rcut;
+
+sph_scan = 0.00001;
+
+npts_plot_theta = 25;
+npts_plot_phi = 50;
+
+p_max = n_max*n_max+2*n_max;
+siwt=siwt-1;
diff --git a/ElectromagneticScattering/scattering_init.py b/ElectromagneticScattering/scattering_init.py
new file mode 100644
index 0000000000000000000000000000000000000000..2f7621c6a6b329b909510f8b92f14a5f33529dcf
--- /dev/null
+++ b/ElectromagneticScattering/scattering_init.py
@@ -0,0 +1,258 @@
+# -*- coding: utf-8 -*-
+
+"""
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering_init.py //
+///////////////////////////////
+"""
+
+import sys,os
+import numpy as np
+from scipy.special import jv, yv, hankel1, hankel2
+sys.path.append(os.getcwd())
+from matplotlib import cm
+import pylab as pl
+from scattering_tmp import *
+np.set_printoptions(precision=2)
+pi = np.pi
+
+ps = np.linspace(1,p_max,p_max)
+r_sphs = np.linspace(r_sph_min,r_sph_max,nb_cuts)
+k_Out = 2*pi*np.sqrt(epsr_Out_re)/lambda0
+epsr_In = epsr_In_re+1j*epsr_In_im
+epsr_Out = epsr_Out_re+1j*epsr_Out_im
+
+phi_range = np.linspace(sph_scan,2*pi-sph_scan,npts_phi)
+theta_range = np.linspace(sph_scan, pi-sph_scan,npts_theta)
+[phi_sph,theta_sph]=np.meshgrid(phi_range,theta_range)
+cos_theta = np.cos(theta_range)
+sin_theta = np.sin(phi_range)
+
+def field_VSH_expansion(post_filename):
+ m_max = n_max
+ p_max = n_max**2 +2*n_max
+ FF_Xnm_t = np.zeros((npts_theta,npts_phi,p_max),dtype=complex)
+ FF_Xnm_p = np.zeros((npts_theta,npts_phi,p_max),dtype=complex)
+ FF_erCrossXnm_t = np.zeros((npts_theta,npts_phi,p_max),dtype=complex)
+ FF_erCrossXnm_p = np.zeros((npts_theta,npts_phi,p_max),dtype=complex)
+ #####################################
+ ##### sn, pn ,un
+ ##### Brian's recurrence relations
+ B_Pnpms = np.zeros((npts_theta,m_max+1,n_max+1))
+ B_unpms = np.zeros((npts_theta,m_max+1,n_max+1))
+ B_snpms = np.zeros((npts_theta,m_max+1,n_max+1))
+ ### Init
+ B_Pnpms[:,0,0] = np.sqrt(1./(4.*pi))
+ B_unpms[:,0,0] = 0.
+ B_unpms[:,1,1] = -0.25*np.sqrt(3./pi)
+ for k in range(npts_theta):
+ u=cos_theta[k]
+ for n in range(1,n_max+1):
+ B_Pnpms[k,n,n] = -np.sqrt((2.*float(n)+1.)/(2.*float(n))) * np.sqrt(1.-u**2) * B_Pnpms[k,n-1,n-1]
+ B_Pnpms[k,n-1,n] = u*np.sqrt(2.*float(n)+1.) * B_Pnpms[k,n-1,n-1]
+ for n in range(2,n_max+1):
+ B_unpms[k,n,n] = -np.sqrt( (float(n)*(2.*float(n)+1.)) / (2.*(float(n)+1.)*(float(n)-1.)) ) * np.sqrt(1.-u**2) * B_unpms[k,n-1,n-1]
+ B_unpms[k,n-1,n] = np.sqrt( (2.*float(n)+1.)*(float(n)-1.)/(float(n)+1.) ) * u * B_unpms[k,n-1,n-1]
+ for n in range(2,n_max+1):
+ for m in range(n-2+1):
+ B_Pnpms[k,m,n] = np.sqrt((4.*(float(n))**2-1.)/((float(n))**2-(float(m))**2)) * u * B_Pnpms[k,m,n-1] \
+ - np.sqrt( ( (2.*float(n)+1.)*((float(n)-1.)**2-(float(m))**2) ) \
+ / ( (2.*float(n)-3.)*((float(n))**2-(float(m))**2) ) )*B_Pnpms[k,m,n-2]
+
+ B_unpms[k,m,n] = np.sqrt( ((4.*(float(n))**2-1.)*(float(n)-1.))/((float(n)**2-float(m)**2)*(float(n)+1.)) ) * u * B_unpms[k,m,n-1] \
+ - np.sqrt( ((2.*float(n)+1.) * (float(n)-1.) * (float(n)-2.) * (float(n-m)-1.) *(float(n+m)-1.)) \
+ / ((2.*float(n)-3.) * (float(n)**2-float(m)**2)*float(n)*(float(n)+1.)))*B_unpms[k,m,n-2]
+ for n in range(0,n_max+1):
+ m=0
+ B_snpms[k,m,n] = 1./float(m+1)*np.sqrt((float(n+m)+1.)*(float(n-m))) * np.sqrt(1.-u**2) *\
+ B_unpms[k,m+1,n] + u*B_unpms[k,m,n]
+ for n in range(1,n_max+1):
+ for m in range(1,n+1):
+ B_snpms[k,m,n] = float(n)/float(m) * u * B_unpms[k,m,n] - float(m+n)/float(m) * \
+ np.sqrt( ( (2.*float(n)+1.)*(float(n)-float(m))*(float(n)-1.) ) / \
+ ( (2.*float(n)-1.)*(float(n)+float(m))*(float(n)+1.) ) )*B_unpms[k,m,n-1]
+ B_Pnmms = np.zeros_like(B_Pnpms)
+ B_unmms = np.zeros_like(B_unpms)
+ B_snmms = np.zeros_like(B_snpms)
+ for m in range(m_max+1):
+ B_Pnmms[:,m,:] = (-1.0)**m * B_Pnpms[:,m,:]
+ B_unmms[:,m,:] = (-1.0)**(m+1) * B_unpms[:,m,:]
+ B_snmms[:,m,:] = (-1.0)**(m) * B_snpms[:,m,:]
+ B_Pnmms=B_Pnmms[:,::-1,:]
+ B_unmms=B_unmms[:,::-1,:]
+ B_snmms=B_snmms[:,::-1,:]
+ B_Pnms = np.concatenate((B_Pnmms,B_Pnpms[:,1::,:]),axis=1)
+ B_unms = np.concatenate((B_unmms,B_unpms[:,1::,:]),axis=1)
+ B_snms = np.concatenate((B_snmms,B_snpms[:,1::,:]),axis=1)
+
+ #####################################
+ ##### sn, pn ,un
+ ##### Brian's recurrence relations
+ m_max = n_max
+ p_max = n_max*n_max+2*n_max
+ aM_nm = np.zeros(p_max,dtype=complex)
+ bN_nm = np.zeros(p_max,dtype=complex)
+ fenm_Y = np.zeros(p_max,dtype=complex)
+ fenm_Z = np.zeros(p_max,dtype=complex)
+ fhnm_X = np.zeros(p_max,dtype=complex)
+
+ B_Pnpms = np.zeros((npts_theta,m_max+1,n_max+1))
+ B_unpms = np.zeros((npts_theta,m_max+1,n_max+1))
+ B_snpms = np.zeros((npts_theta,m_max+1,n_max+1))
+ ### Init
+ B_Pnpms[:,0,0] = np.sqrt(1./(4.*pi))
+ B_unpms[:,0,0] = 0.
+ B_unpms[:,1,1] = -0.25*np.sqrt(3./pi)
+ for k in range(npts_theta):
+ u=cos_theta[k]
+ for n in range(1,n_max+1):
+ B_Pnpms[k,n,n] = -np.sqrt((2.*float(n)+1.)/(2.*float(n))) * np.sqrt(1.-u**2) * B_Pnpms[k,n-1,n-1]
+ B_Pnpms[k,n-1,n] = u*np.sqrt(2.*float(n)+1.) * B_Pnpms[k,n-1,n-1]
+ for n in range(2,n_max+1):
+ B_unpms[k,n,n] = -np.sqrt( (float(n)*(2.*float(n)+1.)) / (2.*(float(n)+1.)*(float(n)-1.)) ) * np.sqrt(1.-u**2) * B_unpms[k,n-1,n-1]
+ B_unpms[k,n-1,n] = np.sqrt( (2.*float(n)+1.)*(float(n)-1.)/(float(n)+1.) ) * u * B_unpms[k,n-1,n-1]
+ for n in range(2,n_max+1):
+ for m in range(n-2+1):
+ B_Pnpms[k,m,n] = np.sqrt((4.*(float(n))**2-1.)/((float(n))**2-(float(m))**2)) * u * B_Pnpms[k,m,n-1] \
+ - np.sqrt( ( (2.*float(n)+1.)*((float(n)-1.)**2-(float(m))**2) ) \
+ / ( (2.*float(n)-3.)*((float(n))**2-(float(m))**2) ) )*B_Pnpms[k,m,n-2]
+
+ B_unpms[k,m,n] = np.sqrt( ((4.*(float(n))**2-1.)*(float(n)-1.))/((float(n)**2-float(m)**2)*(float(n)+1.)) ) * u * B_unpms[k,m,n-1] \
+ - np.sqrt( ((2.*float(n)+1.) * (float(n)-1.) * (float(n)-2.) * (float(n-m)-1.) *(float(n+m)-1.)) \
+ / ((2.*float(n)-3.) * (float(n)**2-float(m)**2)*float(n)*(float(n)+1.)))*B_unpms[k,m,n-2]
+ for n in range(0,n_max+1):
+ m=0
+ B_snpms[k,m,n] = 1./float(m+1)*np.sqrt((float(n+m)+1.)*(float(n-m))) * np.sqrt(1.-u**2) * B_unpms[k,m+1,n] + u*B_unpms[k,m,n]
+ for n in range(1,n_max+1):
+ for m in range(1,n+1):
+ B_snpms[k,m,n] = float(n)/float(m) * u * B_unpms[k,m,n] - float(m+n)/float(m) * \
+ np.sqrt( ( (2.*float(n)+1.)*(float(n)-float(m))*(float(n)-1.) ) / ( (2.*float(n)-1.)*(float(n)+float(m))*(float(n)+1.) ) )*B_unpms[k,m,n-1]
+ B_Pnmms = np.zeros_like(B_Pnpms)
+ B_unmms = np.zeros_like(B_unpms)
+ B_snmms = np.zeros_like(B_snpms)
+ for m in range(m_max+1):
+ B_Pnmms[:,m,:] = (-1.0)**m * B_Pnpms[:,m,:]
+ B_unmms[:,m,:] = (-1.0)**(m+1) * B_unpms[:,m,:]
+ B_snmms[:,m,:] = (-1.0)**(m) * B_snpms[:,m,:]
+ B_Pnmms=B_Pnmms[:,::-1,:]
+ B_unmms=B_unmms[:,::-1,:]
+ B_snmms=B_snmms[:,::-1,:]
+ B_Pnms = np.concatenate((B_Pnmms,B_Pnpms[:,1::,:]),axis=1)
+ B_unms = np.concatenate((B_unmms,B_unpms[:,1::,:]),axis=1)
+ B_snms = np.concatenate((B_snmms,B_snpms[:,1::,:]),axis=1)
+ E_scat_onsphere_sph = np.array( [np.loadtxt(post_filename,usecols=[8])
+ + 1j*np.loadtxt(post_filename,usecols=[11]),
+ np.loadtxt(post_filename,usecols=[9])
+ + 1j*np.loadtxt(post_filename,usecols=[12]),
+ np.loadtxt(post_filename,usecols=[10])
+ + 1j*np.loadtxt(post_filename,usecols=[13])])
+ E_scat_onsphere_sph_r = E_scat_onsphere_sph[0,:].reshape(npts_phi,npts_theta,order='F').transpose()
+ E_scat_onsphere_sph_t = E_scat_onsphere_sph[1,:].reshape(npts_phi,npts_theta,order='F').transpose()
+ E_scat_onsphere_sph_p = E_scat_onsphere_sph[2,:].reshape(npts_phi,npts_theta,order='F').transpose()
+ for ko in range(p_max):
+ po = ps[ko]
+ n = int(np.sqrt(po))
+ m = n*(n+1) - int(po)
+ # print('=========>> po',po,'n',n,'m',m)
+ B_PnmN_costheta = np.tile( B_Pnms[:,m+m_max,n],(npts_phi,1)).transpose()
+ B_UnmN_costheta = np.tile( B_unms[:,m+m_max,n],(npts_phi,1)).transpose()
+ B_SnmN_costheta = np.tile( B_snms[:,m+m_max,n],(npts_phi,1)).transpose()
+ B_Ynm_r = B_PnmN_costheta * np.exp(1j*float(m)*phi_sph)
+ B_Ynm_t = np.zeros_like(phi_sph)
+ B_Ynm_p = np.zeros_like(phi_sph)
+ B_Xnm_r = np.zeros_like(phi_sph)
+ B_Xnm_t = 1j * B_UnmN_costheta * np.exp(1j*float(m)*phi_sph)
+ B_Xnm_p = -1. * B_SnmN_costheta * np.exp(1j*float(m)*phi_sph)
+ B_Znm_r = np.zeros_like(phi_sph)
+ B_Znm_t = 1. * B_SnmN_costheta * np.exp(1j*float(m)*phi_sph)
+ B_Znm_p = 1j * B_UnmN_costheta * np.exp(1j*float(m)*phi_sph)
+ B_erCrossXnm_r = np.zeros_like(phi_sph,dtype=complex)
+ B_erCrossXnm_t = -B_Xnm_p
+ B_erCrossXnm_p = B_Xnm_t
+ FF_Xnm_t[:,:,ko] = B_Xnm_t
+ FF_Xnm_p[:,:,ko] = B_Xnm_p
+ FF_erCrossXnm_t[:,:,ko] = B_erCrossXnm_t
+ FF_erCrossXnm_p[:,:,ko] = B_erCrossXnm_p
+
+
+ sph_bessel_n_ofkr = np.sqrt(pi/(2.*k_Out*r_sph))*outgoing_sph_hankel(float(n )+0.5,k_Out*r_sph)
+ sph_bessel_nminus1_ofkr = np.sqrt(pi/(2.*k_Out*r_sph))*outgoing_sph_hankel(float(n-1)+0.5,k_Out*r_sph)
+ dRicatti_dx_ofkr = (k_Out * r_sph * (sph_bessel_nminus1_ofkr-(float(n+1)/((k_Out*r_sph))) * sph_bessel_n_ofkr) + sph_bessel_n_ofkr)
+ B_Mnm_r = 0.
+ B_Mnm_t = sph_bessel_n_ofkr * B_Xnm_t
+ B_Mnm_p = sph_bessel_n_ofkr * B_Xnm_p
+ B_Nnm_r = 1./(k_Out*r_sph) * np.sqrt(float(n*(n+1))) * sph_bessel_n_ofkr * B_Ynm_r
+ B_Nnm_t = 1./(k_Out*r_sph) * dRicatti_dx_ofkr * B_Znm_t
+ B_Nnm_p = 1./(k_Out*r_sph) * dRicatti_dx_ofkr * B_Znm_p
+
+ B_EdotconjYnm = E_scat_onsphere_sph_r*B_Ynm_r.conjugate()
+ B_EdotconjZnm = E_scat_onsphere_sph_t*B_Znm_t.conjugate() + E_scat_onsphere_sph_p*B_Znm_p.conjugate()
+ B_EdotconjXnm = E_scat_onsphere_sph_t*B_Xnm_t.conjugate() + E_scat_onsphere_sph_p*B_Xnm_p.conjugate()
+ normalize_fhnm_X = 1./sph_bessel_n_ofkr
+ normalize_fenm_Y = k_Out*r_sph/(sph_bessel_n_ofkr*np.sqrt(float(n)*(float(n)+1.)) )
+ normalize_fenm_Z = k_Out*r_sph/dRicatti_dx_ofkr
+ fenm_Y[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjYnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fenm_Y
+ fenm_Z[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjZnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fenm_Z
+ fhnm_X[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fhnm_X
+
+ EdotconjMnm = E_scat_onsphere_sph_r*B_Mnm_r.conjugate() + E_scat_onsphere_sph_t*B_Mnm_t.conjugate() + E_scat_onsphere_sph_p*B_Mnm_p.conjugate()
+ EdotconjNnm = E_scat_onsphere_sph_r*B_Nnm_r.conjugate() + E_scat_onsphere_sph_t*B_Nnm_t.conjugate() + E_scat_onsphere_sph_p*B_Nnm_p.conjugate()
+ MnmconjMnm = B_Mnm_r*B_Mnm_r.conjugate() + B_Mnm_t*B_Mnm_t.conjugate() + B_Mnm_p*B_Mnm_p.conjugate()
+ NnmconjNnm = B_Nnm_r*B_Nnm_r.conjugate() + B_Nnm_t*B_Nnm_t.conjugate() + B_Nnm_p*B_Nnm_p.conjugate()
+ MnmconjNnm = B_Mnm_r*B_Nnm_r.conjugate() + B_Mnm_t*B_Nnm_t.conjugate() + B_Mnm_p*B_Nnm_p.conjugate()
+ XnmconjXnm = B_Xnm_r*B_Xnm_r.conjugate() + B_Xnm_t*B_Xnm_t.conjugate() + B_Xnm_p*B_Xnm_p.conjugate()
+ YnmconjYnm = B_Ynm_r*B_Ynm_r.conjugate() + B_Ynm_t*B_Ynm_t.conjugate() + B_Ynm_p*B_Ynm_p.conjugate()
+ ZnmconjZnm = B_Znm_r*B_Znm_r.conjugate() + B_Znm_t*B_Znm_t.conjugate() + B_Znm_p*B_Znm_p.conjugate()
+ XnmconjYnm = B_Xnm_r*B_Ynm_r.conjugate() + B_Xnm_t*B_Ynm_t.conjugate() + B_Xnm_p*B_Ynm_p.conjugate()
+ YnmconjZnm = B_Ynm_r*B_Znm_r.conjugate() + B_Ynm_t*B_Znm_t.conjugate() + B_Ynm_p*B_Znm_p.conjugate()
+ ZnmconjXnm = B_Znm_r*B_Xnm_r.conjugate() + B_Znm_t*B_Xnm_t.conjugate() + B_Znm_p*B_Xnm_p.conjugate()
+ normalize_aM_nm2 = np.trapz(np.trapz((np.sin(theta_sph)*MnmconjMnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ normalize_bN_nm2 = np.trapz(np.trapz((np.sin(theta_sph)*NnmconjNnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orth = np.trapz(np.trapz((np.sin(theta_sph)*MnmconjNnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orthX = np.trapz(np.trapz((np.sin(theta_sph)*XnmconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orthY = np.trapz(np.trapz((np.sin(theta_sph)*YnmconjYnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orthZ = np.trapz(np.trapz((np.sin(theta_sph)*ZnmconjZnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orth1 = np.trapz(np.trapz((np.sin(theta_sph)*XnmconjYnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orth2 = np.trapz(np.trapz((np.sin(theta_sph)*YnmconjZnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orth3 = np.trapz(np.trapz((np.sin(theta_sph)*ZnmconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ aM_nm[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*EdotconjMnm).transpose(),theta_sph[:,0]),phi_sph[0,:]) / normalize_aM_nm2
+ bN_nm[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*EdotconjNnm).transpose(),theta_sph[:,0]),phi_sph[0,:]) / normalize_bN_nm2
+ return [fhnm_X,fenm_Y,fenm_Z,aM_nm,bN_nm,FF_Xnm_t,FF_Xnm_p,FF_erCrossXnm_t,FF_erCrossXnm_p]
+
+###########################
+def plot_farfield(far_field_sph,filename):
+ vmax_sph = np.max(far_field_sph)
+ vmin_sph = np.min(far_field_sph)
+ x_sph = far_field_sph/vmax_sph*np.sin(theta_sph) * np.cos(phi_sph)
+ y_sph = far_field_sph/vmax_sph*np.sin(theta_sph) * np.sin(phi_sph)
+ z_sph = far_field_sph/vmax_sph*np.cos(theta_sph)
+ fig = pl.figure()
+ ax = fig.add_subplot(111, projection='3d')
+ ax.set_aspect('equal')
+ surf=ax.plot_surface(x_sph,y_sph,z_sph,\
+ facecolors=cm.viridis(far_field_sph/vmax_sph),\
+ rstride=2,\
+ cstride=2,\
+ linewidth=0,\
+ vmin = vmin_sph,\
+ vmax = vmax_sph,\
+ shade=True,\
+ alpha=0.5,\
+ antialiased=False)
+ cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='z', offset=-1)
+ cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='x', offset=-1)
+ cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='y', offset=1 )
+ surf.set_edgecolor('k')
+ max_range = 0.5*np.max(np.array([x_sph.max()-x_sph.min(), y_sph.max()-y_sph.min(), z_sph.max()-z_sph.min()]))
+ mid_x = (x_sph.max()+x_sph.min())*0.5
+ mid_y = (y_sph.max()+y_sph.min())*0.5
+ mid_z = (z_sph.max()+z_sph.min())*0.5
+ ax.set_xlim(mid_x-max_range, mid_x+max_range)
+ ax.set_ylim(mid_y-max_range, mid_y+max_range)
+ ax.set_zlim(mid_z-max_range, mid_z+max_range)
+ ax.xaxis.set_ticklabels([]);ax.xaxis.set_label('x')
+ ax.yaxis.set_ticklabels([]);ax.yaxis.set_label('y')
+ ax.zaxis.set_ticklabels([]);ax.zaxis.set_label('z')
+ pl.savefig(filename,bbox_inches='tight')
+ pl.close('all')
diff --git a/ElectromagneticScattering/scattering_post.py b/ElectromagneticScattering/scattering_post.py
new file mode 100644
index 0000000000000000000000000000000000000000..7e52c8538578872efa6ff9abdab109813e6ef708
--- /dev/null
+++ b/ElectromagneticScattering/scattering_post.py
@@ -0,0 +1,166 @@
+"""
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering_post.py //
+///////////////////////////////
+"""
+
+Tmatrix_rfull_ab = np.zeros((2*p_max,2*p_max,nb_cuts),dtype=complex)
+Tmatrix_rfull_fy = np.zeros((2*p_max,2*p_max,nb_cuts),dtype=complex)
+Tmatrix_rfull_fz = np.zeros((2*p_max,2*p_max,nb_cuts),dtype=complex)
+tab_E_NTF_t_G = np.zeros((npts_theta,npts_phi,3),dtype=complex)
+tab_E_NTF_p_G = np.zeros((npts_theta,npts_phi,3),dtype=complex)
+tab_E_NTF_t_PW = np.zeros((npts_theta,npts_phi),dtype=complex)
+tab_E_NTF_p_PW = np.zeros((npts_theta,npts_phi),dtype=complex)
+
+my_dir='./run_results/'
+if flag_study==0:
+ nbNM = 1
+ p_max_in=1
+elif flag_study==1:
+ nbNM = 2
+ p_max_in=p_max
+elif flag_study==2:
+ nbNM = 1
+ p_max_in=p_max
+ nb_cuts = 3
+elif flag_study==3:
+ p_max_in=0
+
+if siwt==1:
+ outgoing_sph_hankel = hankel2
+else:
+ outgoing_sph_hankel = hankel1
+
+print('Postprocessing...')
+for ke in range(p_max_in):
+ for isN in range(nbNM):
+ pe = ps[ke]
+ ne =int(np.sqrt(pe))
+ me = ne*(ne+1) - int(pe)
+ aM_nm = np.zeros((nb_cuts,p_max),dtype=complex)
+ bN_nm = np.zeros((nb_cuts,p_max),dtype=complex)
+ fenm_Y = np.zeros((nb_cuts,p_max),dtype=complex)
+ fenm_Z = np.zeros((nb_cuts,p_max),dtype=complex)
+ fhnm_X = np.zeros((nb_cuts,p_max),dtype=complex)
+ for nr in range(nb_cuts):
+ r_sph = r_sphs[nr]
+ # print('======>> postprocessing r_sph',r_sph)
+ par_gmsh_getdp=open('parameters_gmsh_getdp.dat', 'a')
+ par_gmsh_getdp.write('r_sph = %3.15e;\n' %(r_sph) )
+ par_gmsh_getdp.close()
+ if flag_study==1:
+ if isN==0:post_filename = my_dir+'E_scat_onsphere_sph_M%g_r%g.dat'%((pe,nr))
+ else: post_filename = my_dir+'E_scat_onsphere_sph_N%g_r%g.dat'%((pe,nr))
+ elif flag_study==0:
+ post_filename = my_dir+'E_scat_onsphere_sph_PW_r%g.dat'%(nr)
+ elif flag_study==2:
+ post_filename = my_dir+'E_tot_onsphere_sph_G%g.dat'%(nr)
+ res=field_VSH_expansion(post_filename)
+ fhnm_X[nr,:] = res[0]
+ fenm_Y[nr,:] = res[1]
+ fenm_Z[nr,:] = res[2]
+ aM_nm[nr,:] = res[3]
+ bN_nm[nr,:] = res[4]
+ FF_Xnm_t = res[5]
+ FF_Xnm_p = res[6]
+ FF_erCrossXnm_t = res[7]
+ FF_erCrossXnm_p = res[8]
+ if flag_study==1:
+ if isN==0:
+ Tmatrix_rfull_ab[ke,0:p_max,:] = aM_nm.transpose()
+ Tmatrix_rfull_ab[ke,p_max:2*p_max,:] = bN_nm.transpose()
+ Tmatrix_rfull_fy[ke,0:p_max,:] = fhnm_X.transpose()
+ Tmatrix_rfull_fy[ke,p_max:2*p_max,:] = fenm_Y.transpose()
+ Tmatrix_rfull_fz[ke,0:p_max,:] = fhnm_X.transpose()
+ Tmatrix_rfull_fz[ke,p_max:2*p_max,:] = fenm_Z.transpose()
+ if isN==1:
+ Tmatrix_rfull_ab[p_max+ke,0:p_max,:] = aM_nm.transpose()
+ Tmatrix_rfull_ab[p_max+ke,p_max:2*p_max,:] = bN_nm.transpose()
+ Tmatrix_rfull_fy[p_max+ke,0:p_max,:] = fhnm_X.transpose()
+ Tmatrix_rfull_fy[p_max+ke,p_max:2*p_max,:] = fenm_Y.transpose()
+ Tmatrix_rfull_fz[p_max+ke,0:p_max,:] = fhnm_X.transpose()
+ Tmatrix_rfull_fz[p_max+ke,p_max:2*p_max,:] = fenm_Z.transpose()
+ Tmatrix_ab = np.mean(Tmatrix_rfull_ab,axis=2)
+ Tmatrix_fy = np.mean(Tmatrix_rfull_fy,axis=2)
+ Tmatrix_fz = np.mean(Tmatrix_rfull_fz,axis=2)
+ std_Tmatrix_ab = np.std(Tmatrix_rfull_ab,axis=2)
+ std_Tmatrix_fy = np.std(Tmatrix_rfull_fy,axis=2)
+ std_Tmatrix_fz = np.std(Tmatrix_rfull_fz,axis=2)
+ np.savez('T_matrix_ellips_sphPML_epsre_%.2lf_eps_im_%.2lf.npz'%(epsr_In.real,epsr_In.imag), \
+ Tmatrix_rfull_ab = Tmatrix_rfull_ab, \
+ Tmatrix_rfull_fy = Tmatrix_rfull_fy, \
+ Tmatrix_rfull_fz = Tmatrix_rfull_fz,\
+ r_sphs = r_sphs,\
+ n_max = n_max,\
+ ps = ps,\
+ p_max = p_max,\
+ lambda0 = lambda0, \
+ Tmatrix_ab = Tmatrix_ab, \
+ Tmatrix_fy = Tmatrix_fy, \
+ Tmatrix_fz = Tmatrix_fz, \
+ std_Tmatrix_ab = std_Tmatrix_ab,\
+ std_Tmatrix_fy = std_Tmatrix_fy,\
+ std_Tmatrix_fz = std_Tmatrix_fz)
+ if flag_study==0:
+ mean_aM_nm = np.mean(aM_nm ,axis=0)
+ mean_bN_nm = np.mean(bN_nm ,axis=0)
+ mean_fenm_Y = np.mean(fenm_Y,axis=0)
+ mean_fenm_Z = np.mean(fenm_Z,axis=0)
+ mean_fhnm_X = np.mean(fhnm_X,axis=0)
+ std_aM_nm = np.std( aM_nm ,axis=0)
+ std_bN_nm = np.std( bN_nm ,axis=0)
+ std_fenm_Y = np.std( fenm_Y,axis=0)
+ std_fenm_Z = np.std( fenm_Z,axis=0)
+ std_fhnm_X = np.std( fhnm_X,axis=0)
+ for ko in range(p_max):
+ tab_E_NTF_t_PW += (1j)**ko * mean_aM_nm[ko]*1j*FF_Xnm_t[:,:,ko] + (1j)**ko * mean_bN_nm[ko]*FF_erCrossXnm_t[:,:,ko]
+ tab_E_NTF_p_PW += (1j)**ko * mean_aM_nm[ko]*1j*FF_Xnm_p[:,:,ko] + (1j)**ko * mean_bN_nm[ko]*FF_erCrossXnm_p[:,:,ko]
+ I_NTF_PW = (np.abs(tab_E_NTF_t_PW))**2 + (np.abs(tab_E_NTF_p_PW))**2
+ plot_farfield(I_NTF_PW[:,:],my_dir+'farfield_PW.pdf')
+ if flag_study==2:
+ tens_Green_scat=np.zeros((3,3))
+ for nr in range(nb_cuts):
+ tens_Green_scat[:,nr] = np.loadtxt(my_dir+'E_scat_im_onpt_G%g.dat'%(nr))[8:11]
+ for ko in range(p_max):
+ tab_E_NTF_t_G[:,:,nr] += (1j)**ko * aM_nm[nr,ko]*1j*FF_Xnm_t[:,:,ko] + (1j)**ko * bN_nm[nr,ko]*FF_erCrossXnm_t[:,:,ko]
+ tab_E_NTF_p_G[:,:,nr] += (1j)**ko * aM_nm[nr,ko]*1j*FF_Xnm_p[:,:,ko] + (1j)**ko * bN_nm[nr,ko]*FF_erCrossXnm_p[:,:,ko]
+ I_NTF_G = (np.abs(tab_E_NTF_t_G))**2 + (np.abs(tab_E_NTF_p_G))**2
+ tens_Green_inc = -siwt*np.eye(3)*k_Out/(6.*pi)
+ tens_Green_tot = tens_Green_scat+tens_Green_inc
+ for k in range(3):plot_farfield(I_NTF_G[:,:,k],my_dir+'farfield_green%g.pdf'%(k))
+
+if flag_study==1:
+ print('Tmatrix_ab:')
+ print(Tmatrix_ab)
+ print('std_Tmatrix_ab:')
+ print(std_Tmatrix_ab)
+ print('Tmatrix_fy:')
+ print(Tmatrix_fy)
+ print('np.abs(2*Tmatrix_fy+1):')
+ print(np.abs(2*Tmatrix_fy+1))
+ print('Tmatrix_fz:')
+ print(Tmatrix_fz)
+ print('np.abs(2*Tmatrix_fz+1):')
+ print(np.abs(2*Tmatrix_fz+1))
+ print('np.abs(2*Tmatrix_ab+1):')
+ print(np.abs(2*Tmatrix_ab+1))
+elif flag_study==0:
+ print('anm:')
+ print(mean_aM_nm)
+ print('bnm')
+ print(mean_bN_nm)
+elif flag_study==2:
+ print('tens_Green_tot')
+ print(tens_Green_tot)
+ print('tens_Green_tot/G0')
+ print(tens_Green_tot/tens_Green_inc[0,0])
+
+elif flag_study==3:
+ eigs = np.loadtxt(my_dir+'EigenValuesReal.txt',usecols=[5]) + 1j*np.loadtxt(my_dir+'EigenValuesImag.txt',usecols=[5])
+ eigs *= (nm/1e-9)**2
+ pl.savez(my_dir+'eigs.npz',eigs=eigs)
+ pl.figure(figsize=(20,15))
+ pl.plot( eigs.real , eigs.imag,'+',mfc='none')
+ pl.grid()
+ pl.savefig(my_dir+'eigenvalues.pdf')
\ No newline at end of file
diff --git a/ElectromagneticScattering/screenshot1_512.png b/ElectromagneticScattering/screenshot1_512.png
new file mode 100644
index 0000000000000000000000000000000000000000..410e93a1a5a34c81866c3761356dd145698aeafe
Binary files /dev/null and b/ElectromagneticScattering/screenshot1_512.png differ