Skip to content
Snippets Groups Projects
Commit 51a10db4 authored by Guillaume Demesy's avatar Guillaume Demesy
Browse files

Adding ElectromagneticScattering model

parent 5031e7ee
No related branches found
No related tags found
No related merge requests found
Pipeline #
Show whitespace changes
Inline Side-by-side
Showing
with 2091 additions and 0 deletions
ElectromagneticScattering/README.md 0 → 100644
+ 23
− 0
View file @ 51a10db4
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
ElectromagneticScattering/scattering.geo 0 → 100644
+ 753
− 0
View file @ 51a10db4
///////////////////////////////
// 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";
ElectromagneticScattering/scattering.pro 0 → 100644
+ 659
− 0
View file @ 51a10db4
///////////////////////////////
// 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
ElectromagneticScattering/scattering_data.geo 0 → 100644
+ 215
− 0
View file @ 51a10db4
///////////////////////////////
// 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;
ElectromagneticScattering/scattering_init.py 0 → 100644
+ 258
− 0
View file @ 51a10db4
# -*- 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')
ElectromagneticScattering/scattering_post.py 0 → 100644
+ 166
− 0
View file @ 51a10db4
"""
///////////////////////////////
// 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
ElectromagneticScattering/scattering_tmp.py 0 → 100644
+ 17
− 0
View file @ 51a10db4
nm = 1.000000000e+00;
lambda0 = 5.876000000e+02;
rbb = 7.345000000e+01;
epsr_In_re = 9.000000000e+00;
epsr_In_im = 0.000000000e+00;
epsr_Out_re = 1.000000000e+00;
epsr_Out_im = 0.000000000e+00;
sph_scan = 1.000000000e-05;
r_sph_min = 1.322100000e+02;
r_sph_max = 2.497300000e+02;
nb_cuts = 3;
npts_theta = 50;
npts_phi = 100;
flag_study = 2;
n_max = 1;
p_max = 3;
siwt = -1;
ElectromagneticScattering/screenshot1_512.png 0 → 100644
+ 0
− 0
View file @ 51a10db4
ElectromagneticScattering/screenshot1_512.png

895 KiB

0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment