Main Lemma Repository
選択できるのは25トピックまでです。 トピックは、先頭が英数字で、英数字とダッシュ('-')を使用した35文字以内のものにしてください。

DipoleSource.cpp 74KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334
  1. /* This file is part of Lemma, a geophysical modelling and inversion API */
  2. /* This Source Code Form is subject to the terms of the Mozilla Public
  3. * License, v. 2.0. If a copy of the MPL was not distributed with this
  4. * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
  5. /**
  6. @file
  7. @author Trevor Irons
  8. @date 12/02/2009
  9. **/
  10. #include "DipoleSource.h"
  11. #include "KernelEM1DManager.h"
  12. //#include "GroundedElectricDipole.h"
  13. //#include "UngroundedElectricDipole.h"
  14. //#include "MagneticDipole.h"
  15. #include "FieldPoints.h"
  16. #include "HankelTransform.h"
  17. namespace Lemma {
  18. // ==================== FRIENDS ======================
  19. std::ostream &operator<<(std::ostream &stream, const DipoleSource &ob) {
  20. stream << ob.Serialize() << "\n";
  21. return stream;
  22. }
  23. /*
  24. bool DipoleSource::operator == (DipoleSource& rhs) const {
  25. if (Location != rhs.Location) return false;
  26. return true;
  27. }
  28. */
  29. // ==================== LIFECYCLE ======================
  30. DipoleSource::DipoleSource( const ctor_key& key ) : LemmaObject( key ),
  31. Type(NOSOURCETYPE),
  32. irec(-1),
  33. Phase(0),
  34. Moment(1),
  35. KernelManager(nullptr),
  36. Receivers(nullptr),
  37. Earth(nullptr)
  38. {
  39. this->Location.setZero();
  40. this->Phat.setZero();
  41. }
  42. DipoleSource::DipoleSource( const YAML::Node& node, const ctor_key& key ) : LemmaObject( node, key ),
  43. Type(NOSOURCETYPE),
  44. irec(-1),
  45. Phase(0),
  46. Moment(1),
  47. KernelManager(nullptr),
  48. Receivers(nullptr),
  49. Earth(nullptr)
  50. {
  51. Type = string2Enum<DIPOLESOURCETYPE>(node["Type"].as<std::string>());
  52. this->Location = node["Location"].as<Vector3r>();
  53. this->Phat.setZero();
  54. }
  55. DipoleSource::~DipoleSource() {
  56. }
  57. std::shared_ptr<DipoleSource> DipoleSource::NewSP() {
  58. return std::make_shared<DipoleSource> ( ctor_key() );
  59. }
  60. YAML::Node DipoleSource::Serialize() const {
  61. YAML::Node node = LemmaObject::Serialize();
  62. node.SetTag( GetName() );
  63. node["Type"] = enum2String(Type);
  64. node["Location"] = Location;
  65. node["Phat"] = Phat;
  66. node["Freqs"] = Freqs;
  67. node["Phase"] = Phase;
  68. node["Moment"] = Moment;
  69. return node;
  70. }
  71. std::shared_ptr< DipoleSource > DipoleSource::DeSerialize(const YAML::Node& node) {
  72. if (node.Tag() != "DipoleSource") {
  73. throw DeSerializeTypeMismatch( "DipoleSource", node.Tag());
  74. }
  75. return std::make_shared<DipoleSource> ( node, ctor_key() );
  76. }
  77. std::shared_ptr<DipoleSource> DipoleSource::Clone() {
  78. auto Obj = DipoleSource::NewSP();
  79. // copy
  80. Obj->Type = Type;
  81. Obj->irec = irec;
  82. Obj->lays = lays;
  83. Obj->layr = layr;
  84. Obj->Phase = Phase;
  85. Obj->Moment = Moment;
  86. Obj->xxp = xxp;
  87. Obj->yyp = yyp;
  88. Obj->rho = rho;
  89. Obj->sp = sp;
  90. Obj->cp = cp;
  91. Obj->scp = scp;
  92. Obj->sps = sps;
  93. Obj->cps = cps;
  94. Obj->c2p = c2p;
  95. Obj->FieldsToCalculate = FieldsToCalculate;
  96. Obj->f = f;
  97. Obj->ik = ik;
  98. Obj->Location = Location;
  99. Obj->Phat = Phat;
  100. Obj->Freqs = Freqs;
  101. return Obj;
  102. }
  103. //--------------------------------------------------------------------------------------
  104. // Class: DipoleSource
  105. // Method: GetName
  106. // Description: Class identifier
  107. //--------------------------------------------------------------------------------------
  108. inline std::string DipoleSource::GetName ( ) const {
  109. return CName;
  110. } // ----- end of method DipoleSource::GetName -----
  111. // ==================== ACCESS ======================
  112. void DipoleSource::SetLocation(const Vector3r &posin) {
  113. this->Location = posin;
  114. }
  115. void DipoleSource::SetLocation(const Real &xp, const Real &yp,
  116. const Real &zp) {
  117. this->Location = Vector3r(xp, yp, zp);
  118. }
  119. void DipoleSource::SetPhase(const Real &phase) {
  120. this->Phase = phase;
  121. }
  122. void DipoleSource::SetPolarity(const DipoleSourcePolarity &pol) {
  123. static bool called = false;
  124. if (!called) {
  125. std::cerr << "\n\n=================================================================\n"
  126. << "WARNING: Use of deprecated method DipoleSource::SetPolarity(pol)\n"
  127. << "Use more general SetPolarisation( Vector3r ) or SetPolarisation( x, y, z );\n"
  128. << "This method will be removed in future versions of Lemma"
  129. << "\n=================================================================\n";
  130. called = true;
  131. }
  132. // Polarity = pol;
  133. // switch (Polarity) {
  134. // case POSITIVE:
  135. // Moment = std::abs(Moment);
  136. // break;
  137. // case NEGATIVE:
  138. // Moment = -std::abs(Moment);
  139. // break;
  140. // default:
  141. // throw NonValidDipolePolarity();
  142. // };
  143. }
  144. void DipoleSource::SetType(const DIPOLESOURCETYPE & stype) {
  145. switch (stype) {
  146. case (GROUNDEDELECTRICDIPOLE):
  147. this->Type = stype;
  148. break;
  149. case (UNGROUNDEDELECTRICDIPOLE):
  150. this->Type = stype;
  151. break;
  152. case (MAGNETICDIPOLE):
  153. this->Type = stype;
  154. break;
  155. default:
  156. throw NonValidDipoleTypeAssignment();
  157. }
  158. }
  159. void DipoleSource::SetPolarisation(const Vector3r& pol) {
  160. this->Phat = pol / pol.norm();
  161. }
  162. void DipoleSource::SetPolarisation(const Real& x, const Real& y, const Real& z) {
  163. Vector3r pol = (VectorXr(3) << x, y, z).finished();
  164. this->Phat = pol / pol.norm();
  165. }
  166. Vector3r DipoleSource::GetPolarisation() {
  167. return Phat;
  168. }
  169. DIPOLESOURCETYPE DipoleSource::GetType() {
  170. return Type;
  171. }
  172. void DipoleSource::SetPolarisation(const
  173. DipoleSourcePolarisation &pol) {
  174. static bool called = false;
  175. if (!called) {
  176. std::cout << "\n\n========================================================================================\n"
  177. << "WARNING: Use of deprecated method DipoleSource::SetPolarisation(DipleSourcePolarisation)\n"
  178. << "Use more general SetPolarisation( Vector3r ) or SetPolarisation( x, y, z );\n"
  179. << "This method will be removed in future versions of Lemma"
  180. << "\n========================================================================================\n";
  181. called = true;
  182. }
  183. switch (pol) {
  184. case (XPOLARISATION):
  185. this->Phat = (VectorXr(3) << 1, 0, 0).finished();
  186. break;
  187. case (YPOLARISATION):
  188. this->Phat = (VectorXr(3) << 0, 1, 0).finished();
  189. break;
  190. case (ZPOLARISATION):
  191. this->Phat = (VectorXr(3) << 0, 0, 1).finished();
  192. break;
  193. default:
  194. throw NonValidDipolePolarisationAssignment();
  195. }
  196. }
  197. void DipoleSource::SetMoment(const Real &moment) {
  198. this->Moment = moment;
  199. }
  200. // ==================== OPERATIONS =====================
  201. void DipoleSource::SetKernels(const int& ifreq, const FIELDCALCULATIONS& Fields , std::shared_ptr<FieldPoints> ReceiversIn, const int& irecin,
  202. std::shared_ptr<LayeredEarthEM> EarthIn ) {
  203. if (Receivers != ReceiversIn) {
  204. Receivers = ReceiversIn;
  205. }
  206. if (Earth != EarthIn) {
  207. Earth = EarthIn;
  208. }
  209. if (irecin != irec) {
  210. irec = irecin;
  211. }
  212. if (FieldsToCalculate != Fields) {
  213. FieldsToCalculate = Fields;
  214. }
  215. xxp = Receivers->GetLocation(irec)[0] - Location[0];
  216. yyp = Receivers->GetLocation(irec)[1] - Location[1];
  217. rho = (Receivers->GetLocation(irec).head<2>() - Location.head<2>()).norm();
  218. sp = yyp/rho;
  219. cp = xxp/rho;
  220. scp = sp*cp;
  221. sps = sp*sp;
  222. cps = cp*cp;
  223. c2p = cps-sps;
  224. f = VectorXcr::Zero(13);
  225. ik = VectorXi::Zero(13);
  226. lays = Earth->GetLayerAtThisDepth(Location[2]);
  227. layr = Earth->GetLayerAtThisDepth(Receivers->GetLocation(irec)[2]);
  228. KernelManager = KernelEM1DManager::NewSP();
  229. KernelManager->SetEarth(Earth);
  230. // alternative is to use weak_ptr here, this is deep and internal, and we are safe.
  231. KernelManager->SetDipoleSource( shared_from_this().get() , ifreq, Receivers->GetLocation(irec)[2]);
  232. kernelFreq = Freqs(ifreq); // this is never used
  233. ReSetKernels( ifreq, Fields, Receivers, irec, Earth );
  234. return;
  235. }
  236. // TODO we could make the dipoles template specializations avoiding this rats nest of switch statements. Probably
  237. // not the most critical piece though
  238. void DipoleSource::ReSetKernels(const int& ifreq, const FIELDCALCULATIONS& Fields , std::shared_ptr<FieldPoints> Receivers, const int& irec,
  239. std::shared_ptr<LayeredEarthEM> Earth ) {
  240. Vector3r Pol = Phat;
  241. switch (Type) {
  242. case (GROUNDEDELECTRICDIPOLE):
  243. if (std::abs(Pol[2]) > 0) { // z dipole
  244. switch(FieldsToCalculate) {
  245. case E:
  246. if (lays == 0 && layr == 0) {
  247. ik[10] = KernelManager->AddKernel<TM, 10, INAIR, INAIR>( );
  248. ik[11] = KernelManager->AddKernel<TM, 11, INAIR, INAIR>( );
  249. } else if (lays == 0 && layr > 0) {
  250. ik[10] = KernelManager->AddKernel<TM, 10, INAIR, INGROUND>( );
  251. ik[11] = KernelManager->AddKernel<TM, 11, INAIR, INGROUND>( );
  252. } else if (lays > 0 && layr == 0) {
  253. ik[10] = KernelManager->AddKernel<TM, 10, INGROUND, INAIR>( );
  254. ik[11] = KernelManager->AddKernel<TM, 11, INGROUND, INAIR>( );
  255. } else {
  256. ik[10] = KernelManager->AddKernel<TM, 10, INGROUND, INGROUND>( );
  257. ik[11] = KernelManager->AddKernel<TM, 11, INGROUND, INGROUND>( );
  258. }
  259. break;
  260. case H:
  261. if (lays == 0 && layr == 0) {
  262. ik[12] = KernelManager->AddKernel<TM, 12, INAIR, INAIR>( );
  263. } else if (lays == 0 && layr > 0) {
  264. ik[12] = KernelManager->AddKernel<TM, 12, INAIR, INGROUND>( );
  265. } else if (lays > 0 && layr == 0) {
  266. ik[12] = KernelManager->AddKernel<TM, 12, INGROUND, INAIR>( );
  267. } else {
  268. ik[12] = KernelManager->AddKernel<TM, 12, INGROUND, INGROUND>( );
  269. }
  270. break;
  271. case BOTH:
  272. if ( lays == 0 && layr == 0) {
  273. ik[10] = KernelManager->AddKernel<TM, 10, INAIR, INAIR>( );
  274. ik[11] = KernelManager->AddKernel<TM, 11, INAIR, INAIR>( );
  275. ik[12] = KernelManager->AddKernel<TM, 12, INAIR, INAIR>( );
  276. } else if (lays == 0 && layr > 0) {
  277. ik[10] = KernelManager->AddKernel<TM, 10, INAIR, INGROUND>( );
  278. ik[11] = KernelManager->AddKernel<TM, 11, INAIR, INGROUND>( );
  279. ik[12] = KernelManager->AddKernel<TM, 12, INAIR, INGROUND>( );
  280. } else if (lays > 0 && layr == 0) {
  281. ik[10] = KernelManager->AddKernel<TM, 10, INGROUND, INAIR>( );
  282. ik[11] = KernelManager->AddKernel<TM, 11, INGROUND, INAIR>( );
  283. ik[12] = KernelManager->AddKernel<TM, 12, INGROUND, INAIR>( );
  284. } else {
  285. ik[10] = KernelManager->AddKernel<TM, 10, INGROUND, INGROUND>( );
  286. ik[11] = KernelManager->AddKernel<TM, 11, INGROUND, INGROUND>( );
  287. ik[12] = KernelManager->AddKernel<TM, 12, INGROUND, INGROUND>( );
  288. }
  289. }
  290. }
  291. if (std::abs(Pol[1]) > 0 || std::abs(Pol[0]) > 0) { // x or y grounded HED dipole
  292. switch(FieldsToCalculate) {
  293. case E:
  294. if ( lays == 0 && layr == 0) {
  295. ik[0] = KernelManager->AddKernel<TM, 0, INAIR, INAIR>( );
  296. ik[1] = KernelManager->AddKernel<TM, 1, INAIR, INAIR>( );
  297. ik[4] = KernelManager->AddKernel<TM, 4, INAIR, INAIR>( );
  298. ik[2] = KernelManager->AddKernel<TE, 2, INAIR, INAIR>( );
  299. ik[3] = KernelManager->AddKernel<TE, 3, INAIR, INAIR>( );
  300. } else if (lays == 0 && layr > 0) {
  301. ik[0] = KernelManager->AddKernel<TM, 0, INAIR, INGROUND>( );
  302. ik[1] = KernelManager->AddKernel<TM, 1, INAIR, INGROUND>( );
  303. ik[4] = KernelManager->AddKernel<TM, 4, INAIR, INGROUND>( );
  304. ik[2] = KernelManager->AddKernel<TE, 2, INAIR, INGROUND>( );
  305. ik[3] = KernelManager->AddKernel<TE, 3, INAIR, INGROUND>( );
  306. } else if (lays > 0 && layr == 0) {
  307. ik[0] = KernelManager->AddKernel<TM, 0, INGROUND, INAIR>( );
  308. ik[1] = KernelManager->AddKernel<TM, 1, INGROUND, INAIR>( );
  309. ik[4] = KernelManager->AddKernel<TM, 4, INGROUND, INAIR>( );
  310. ik[2] = KernelManager->AddKernel<TE, 2, INGROUND, INAIR>( );
  311. ik[3] = KernelManager->AddKernel<TE, 3, INGROUND, INAIR>( );
  312. } else {
  313. ik[0] = KernelManager->AddKernel<TM, 0, INGROUND, INGROUND>( );
  314. ik[1] = KernelManager->AddKernel<TM, 1, INGROUND, INGROUND>( );
  315. ik[4] = KernelManager->AddKernel<TM, 4, INGROUND, INGROUND>( );
  316. ik[2] = KernelManager->AddKernel<TE, 2, INGROUND, INGROUND>( );
  317. ik[3] = KernelManager->AddKernel<TE, 3, INGROUND, INGROUND>( );
  318. }
  319. break;
  320. case H:
  321. if (lays == 0 && layr == 0) {
  322. ik[5] = KernelManager->AddKernel<TM, 5, INAIR, INAIR>( );
  323. ik[6] = KernelManager->AddKernel<TM, 6, INAIR, INAIR>( );
  324. ik[7] = KernelManager->AddKernel<TE, 7, INAIR, INAIR>( );
  325. ik[8] = KernelManager->AddKernel<TE, 8, INAIR, INAIR>( );
  326. ik[9] = KernelManager->AddKernel<TE, 9, INAIR, INAIR>( );
  327. } else if (lays == 0 && layr > 0) {
  328. ik[5] = KernelManager->AddKernel<TM, 5, INAIR, INGROUND>( );
  329. ik[6] = KernelManager->AddKernel<TM, 6, INAIR, INGROUND>( );
  330. ik[7] = KernelManager->AddKernel<TE, 7, INAIR, INGROUND>( );
  331. ik[8] = KernelManager->AddKernel<TE, 8, INAIR, INGROUND>( );
  332. ik[9] = KernelManager->AddKernel<TE, 9, INAIR, INGROUND>( );
  333. } else if (lays > 0 && layr == 0) {
  334. ik[5] = KernelManager->AddKernel<TM, 5, INGROUND, INAIR>( );
  335. ik[6] = KernelManager->AddKernel<TM, 6, INGROUND, INAIR>( );
  336. ik[7] = KernelManager->AddKernel<TE, 7, INGROUND, INAIR>( );
  337. ik[8] = KernelManager->AddKernel<TE, 8, INGROUND, INAIR>( );
  338. ik[9] = KernelManager->AddKernel<TE, 9, INGROUND, INAIR>( );
  339. } else {
  340. ik[5] = KernelManager->AddKernel<TM, 5, INGROUND, INGROUND>( );
  341. ik[6] = KernelManager->AddKernel<TM, 6, INGROUND, INGROUND>( );
  342. ik[7] = KernelManager->AddKernel<TE, 7, INGROUND, INGROUND>( );
  343. ik[8] = KernelManager->AddKernel<TE, 8, INGROUND, INGROUND>( );
  344. ik[9] = KernelManager->AddKernel<TE, 9, INGROUND, INGROUND>( );
  345. }
  346. break;
  347. case BOTH:
  348. if (lays == 0 && layr == 0) {
  349. ik[0] = KernelManager->AddKernel<TM, 0, INAIR, INAIR>( );
  350. ik[1] = KernelManager->AddKernel<TM, 1, INAIR, INAIR>( );
  351. ik[4] = KernelManager->AddKernel<TM, 4, INAIR, INAIR>( );
  352. ik[2] = KernelManager->AddKernel<TE, 2, INAIR, INAIR>( );
  353. ik[3] = KernelManager->AddKernel<TE, 3, INAIR, INAIR>( );
  354. ik[5] = KernelManager->AddKernel<TM, 5, INAIR, INAIR>( );
  355. ik[6] = KernelManager->AddKernel<TM, 6, INAIR, INAIR>( );
  356. ik[7] = KernelManager->AddKernel<TE, 7, INAIR, INAIR>( );
  357. ik[8] = KernelManager->AddKernel<TE, 8, INAIR, INAIR>( );
  358. ik[9] = KernelManager->AddKernel<TE, 9, INAIR, INAIR>( );
  359. } else if (lays == 0 && layr > 0) {
  360. ik[0] = KernelManager->AddKernel<TM, 0, INAIR, INGROUND>( );
  361. ik[1] = KernelManager->AddKernel<TM, 1, INAIR, INGROUND>( );
  362. ik[4] = KernelManager->AddKernel<TM, 4, INAIR, INGROUND>( );
  363. ik[2] = KernelManager->AddKernel<TE, 2, INAIR, INGROUND>( );
  364. ik[3] = KernelManager->AddKernel<TE, 3, INAIR, INGROUND>( );
  365. ik[5] = KernelManager->AddKernel<TM, 5, INAIR, INGROUND>( );
  366. ik[6] = KernelManager->AddKernel<TM, 6, INAIR, INGROUND>( );
  367. ik[7] = KernelManager->AddKernel<TE, 7, INAIR, INGROUND>( );
  368. ik[8] = KernelManager->AddKernel<TE, 8, INAIR, INGROUND>( );
  369. ik[9] = KernelManager->AddKernel<TE, 9, INAIR, INGROUND>( );
  370. } else if (lays > 0 && layr == 0) {
  371. ik[0] = KernelManager->AddKernel<TM, 0, INGROUND, INAIR>( );
  372. ik[1] = KernelManager->AddKernel<TM, 1, INGROUND, INAIR>( );
  373. ik[4] = KernelManager->AddKernel<TM, 4, INGROUND, INAIR>( );
  374. ik[2] = KernelManager->AddKernel<TE, 2, INGROUND, INAIR>( );
  375. ik[3] = KernelManager->AddKernel<TE, 3, INGROUND, INAIR>( );
  376. ik[5] = KernelManager->AddKernel<TM, 5, INGROUND, INAIR>( );
  377. ik[6] = KernelManager->AddKernel<TM, 6, INGROUND, INAIR>( );
  378. ik[7] = KernelManager->AddKernel<TE, 7, INGROUND, INAIR>( );
  379. ik[8] = KernelManager->AddKernel<TE, 8, INGROUND, INAIR>( );
  380. ik[9] = KernelManager->AddKernel<TE, 9, INGROUND, INAIR>( );
  381. } else {
  382. ik[0] = KernelManager->AddKernel<TM, 0, INGROUND, INGROUND>( );
  383. ik[1] = KernelManager->AddKernel<TM, 1, INGROUND, INGROUND>( );
  384. ik[4] = KernelManager->AddKernel<TM, 4, INGROUND, INGROUND>( );
  385. ik[2] = KernelManager->AddKernel<TE, 2, INGROUND, INGROUND>( );
  386. ik[3] = KernelManager->AddKernel<TE, 3, INGROUND, INGROUND>( );
  387. ik[5] = KernelManager->AddKernel<TM, 5, INGROUND, INGROUND>( );
  388. ik[6] = KernelManager->AddKernel<TM, 6, INGROUND, INGROUND>( );
  389. ik[7] = KernelManager->AddKernel<TE, 7, INGROUND, INGROUND>( );
  390. ik[8] = KernelManager->AddKernel<TE, 8, INGROUND, INGROUND>( );
  391. ik[9] = KernelManager->AddKernel<TE, 9, INGROUND, INGROUND>( );
  392. }
  393. break;
  394. }
  395. }
  396. break;
  397. case (UNGROUNDEDELECTRICDIPOLE):
  398. if (std::abs(Pol[2]) > 0) { // z dipole
  399. switch(FieldsToCalculate) {
  400. case E:
  401. if (lays == 0 && layr == 0) {
  402. ik[11] = KernelManager->AddKernel<TM, 11, INAIR, INAIR>( );
  403. } else if (lays == 0 && layr > 0) {
  404. ik[11] = KernelManager->AddKernel<TM, 11, INAIR, INGROUND>( );
  405. } else if (lays > 0 && layr == 0) {
  406. ik[11] = KernelManager->AddKernel<TM, 11, INGROUND, INAIR>( );
  407. } else {
  408. ik[11] = KernelManager->AddKernel<TM, 11, INGROUND, INGROUND>( );
  409. }
  410. break;
  411. case H:
  412. if (lays == 0 && layr == 0) {
  413. ik[12] = KernelManager->AddKernel<TM, 12, INAIR, INAIR>( );
  414. } else if (lays == 0 && layr > 0) {
  415. ik[12] = KernelManager->AddKernel<TM, 12, INAIR, INGROUND>( );
  416. } else if (lays > 0 && layr == 0) {
  417. ik[12] = KernelManager->AddKernel<TM, 12, INGROUND, INAIR>( );
  418. } else {
  419. ik[12] = KernelManager->AddKernel<TM, 12, INGROUND, INGROUND>( );
  420. }
  421. break;
  422. case BOTH:
  423. if ( lays == 0 && layr == 0) {
  424. ik[11] = KernelManager->AddKernel<TM, 11, INAIR, INAIR>( );
  425. ik[12] = KernelManager->AddKernel<TM, 12, INAIR, INAIR>( );
  426. } else if (lays == 0 && layr > 0) {
  427. ik[11] = KernelManager->AddKernel<TM, 11, INAIR, INGROUND>( );
  428. ik[12] = KernelManager->AddKernel<TM, 12, INAIR, INGROUND>( );
  429. } else if (lays > 0 && layr == 0) {
  430. ik[11] = KernelManager->AddKernel<TM, 11, INGROUND, INAIR>( );
  431. ik[12] = KernelManager->AddKernel<TM, 12, INGROUND, INAIR>( );
  432. } else {
  433. ik[11] = KernelManager->AddKernel<TM, 11, INGROUND, INGROUND>( );
  434. ik[12] = KernelManager->AddKernel<TM, 12, INGROUND, INGROUND>( );
  435. }
  436. }
  437. }
  438. if (std::abs(Pol[1]) > 0 || std::abs(Pol[0]) > 0) { // x or y grounded HED dipole
  439. switch(FieldsToCalculate) {
  440. case E:
  441. if ( lays == 0 && layr == 0) {
  442. ik[2] = KernelManager->AddKernel<TE, 2, INAIR, INAIR>( );
  443. ik[3] = KernelManager->AddKernel<TE, 3, INAIR, INAIR>( );
  444. } else if (lays == 0 && layr > 0) {
  445. ik[2] = KernelManager->AddKernel<TE, 2, INAIR, INGROUND>( );
  446. ik[3] = KernelManager->AddKernel<TE, 3, INAIR, INGROUND>( );
  447. } else if (lays > 0 && layr == 0) {
  448. ik[2] = KernelManager->AddKernel<TE, 2, INGROUND, INAIR>( );
  449. ik[3] = KernelManager->AddKernel<TE, 3, INGROUND, INAIR>( );
  450. } else {
  451. ik[2] = KernelManager->AddKernel<TE, 2, INGROUND, INGROUND>( );
  452. ik[3] = KernelManager->AddKernel<TE, 3, INGROUND, INGROUND>( );
  453. }
  454. break;
  455. case H:
  456. if (lays == 0 && layr == 0) {
  457. ik[5] = KernelManager->AddKernel<TM, 5, INAIR, INAIR>( );
  458. ik[6] = KernelManager->AddKernel<TM, 6, INAIR, INAIR>( );
  459. ik[7] = KernelManager->AddKernel<TE, 7, INAIR, INAIR>( );
  460. ik[8] = KernelManager->AddKernel<TE, 8, INAIR, INAIR>( );
  461. ik[9] = KernelManager->AddKernel<TE, 9, INAIR, INAIR>( );
  462. } else if (lays == 0 && layr > 0) {
  463. ik[5] = KernelManager->AddKernel<TM, 5, INAIR, INGROUND>( );
  464. ik[6] = KernelManager->AddKernel<TM, 6, INAIR, INGROUND>( );
  465. ik[7] = KernelManager->AddKernel<TE, 7, INAIR, INGROUND>( );
  466. ik[8] = KernelManager->AddKernel<TE, 8, INAIR, INGROUND>( );
  467. ik[9] = KernelManager->AddKernel<TE, 9, INAIR, INGROUND>( );
  468. } else if (lays > 0 && layr == 0) {
  469. ik[5] = KernelManager->AddKernel<TM, 5, INGROUND, INAIR>( );
  470. ik[6] = KernelManager->AddKernel<TM, 6, INGROUND, INAIR>( );
  471. ik[7] = KernelManager->AddKernel<TE, 7, INGROUND, INAIR>( );
  472. ik[8] = KernelManager->AddKernel<TE, 8, INGROUND, INAIR>( );
  473. ik[9] = KernelManager->AddKernel<TE, 9, INGROUND, INAIR>( );
  474. } else {
  475. ik[5] = KernelManager->AddKernel<TM, 5, INGROUND, INGROUND>( );
  476. ik[6] = KernelManager->AddKernel<TM, 6, INGROUND, INGROUND>( );
  477. ik[7] = KernelManager->AddKernel<TE, 7, INGROUND, INGROUND>( );
  478. ik[8] = KernelManager->AddKernel<TE, 8, INGROUND, INGROUND>( );
  479. ik[9] = KernelManager->AddKernel<TE, 9, INGROUND, INGROUND>( );
  480. }
  481. break;
  482. case BOTH:
  483. if (lays == 0 && layr == 0) {
  484. ik[5] = KernelManager->AddKernel<TM, 5, INAIR, INAIR>( );
  485. ik[6] = KernelManager->AddKernel<TM, 6, INAIR, INAIR>( );
  486. ik[2] = KernelManager->AddKernel<TE, 2, INAIR, INAIR>( );
  487. ik[3] = KernelManager->AddKernel<TE, 3, INAIR, INAIR>( );
  488. ik[7] = KernelManager->AddKernel<TE, 7, INAIR, INAIR>( );
  489. ik[8] = KernelManager->AddKernel<TE, 8, INAIR, INAIR>( );
  490. ik[9] = KernelManager->AddKernel<TE, 9, INAIR, INAIR>( );
  491. } else if (lays == 0 && layr > 0) {
  492. ik[5] = KernelManager->AddKernel<TM, 5, INAIR, INGROUND>( );
  493. ik[6] = KernelManager->AddKernel<TM, 6, INAIR, INGROUND>( );
  494. ik[2] = KernelManager->AddKernel<TE, 2, INAIR, INGROUND>( );
  495. ik[3] = KernelManager->AddKernel<TE, 3, INAIR, INGROUND>( );
  496. ik[7] = KernelManager->AddKernel<TE, 7, INAIR, INGROUND>( );
  497. ik[8] = KernelManager->AddKernel<TE, 8, INAIR, INGROUND>( );
  498. ik[9] = KernelManager->AddKernel<TE, 9, INAIR, INGROUND>( );
  499. } else if (lays > 0 && layr == 0) {
  500. ik[5] = KernelManager->AddKernel<TM, 5, INGROUND, INAIR>( );
  501. ik[6] = KernelManager->AddKernel<TM, 6, INGROUND, INAIR>( );
  502. ik[2] = KernelManager->AddKernel<TE, 2, INGROUND, INAIR>( );
  503. ik[3] = KernelManager->AddKernel<TE, 3, INGROUND, INAIR>( );
  504. ik[7] = KernelManager->AddKernel<TE, 7, INGROUND, INAIR>( );
  505. ik[8] = KernelManager->AddKernel<TE, 8, INGROUND, INAIR>( );
  506. ik[9] = KernelManager->AddKernel<TE, 9, INGROUND, INAIR>( );
  507. } else {
  508. ik[5] = KernelManager->AddKernel<TM, 5, INGROUND, INGROUND>( );
  509. ik[6] = KernelManager->AddKernel<TM, 6, INGROUND, INGROUND>( );
  510. ik[2] = KernelManager->AddKernel<TE, 2, INGROUND, INGROUND>( );
  511. ik[3] = KernelManager->AddKernel<TE, 3, INGROUND, INGROUND>( );
  512. ik[7] = KernelManager->AddKernel<TE, 7, INGROUND, INGROUND>( );
  513. ik[8] = KernelManager->AddKernel<TE, 8, INGROUND, INGROUND>( );
  514. ik[9] = KernelManager->AddKernel<TE, 9, INGROUND, INGROUND>( );
  515. }
  516. break;
  517. }
  518. }
  519. break;
  520. case (MAGNETICDIPOLE):
  521. if (std::abs(Pol[2]) > 0) { // z dipole
  522. switch (FieldsToCalculate) {
  523. case E:
  524. if (lays == 0 && layr == 0) {
  525. ik[12] = KernelManager->AddKernel<TE, 12, INAIR, INAIR>( );
  526. } else if (lays == 0 && layr > 0) {
  527. ik[12] = KernelManager->AddKernel<TE, 12, INAIR, INGROUND>( );
  528. } else if (lays > 0 && layr == 0) {
  529. ik[12] = KernelManager->AddKernel<TE, 12, INGROUND, INAIR>( );
  530. } else {
  531. ik[12] = KernelManager->AddKernel<TE, 12, INGROUND, INGROUND>( );
  532. }
  533. break;
  534. case H:
  535. if (lays == 0 && layr == 0) {
  536. ik[10] = KernelManager->AddKernel<TE, 10, INAIR, INAIR>( );
  537. ik[11] = KernelManager->AddKernel<TE, 11, INAIR, INAIR>( );
  538. } else if (lays == 0 && layr > 0) {
  539. ik[10] = KernelManager->AddKernel<TE, 10, INAIR, INGROUND>( );
  540. ik[11] = KernelManager->AddKernel<TE, 11, INAIR, INGROUND>( );
  541. } else if (lays > 0 && layr == 0) {
  542. ik[10] = KernelManager->AddKernel<TE, 10, INGROUND, INAIR>( );
  543. ik[11] = KernelManager->AddKernel<TE, 11, INGROUND, INAIR>( );
  544. } else {
  545. ik[10] = KernelManager->AddKernel<TE, 10, INGROUND, INGROUND>( );
  546. ik[11] = KernelManager->AddKernel<TE, 11, INGROUND, INGROUND>( );
  547. }
  548. break;
  549. case BOTH:
  550. if (lays == 0 && layr == 0) {
  551. ik[12] = KernelManager->AddKernel<TE, 12, INAIR, INAIR>( );
  552. ik[10] = KernelManager->AddKernel<TE, 10, INAIR, INAIR>( );
  553. ik[11] = KernelManager->AddKernel<TE, 11, INAIR, INAIR>( );
  554. } else if (lays == 0 && layr > 0) {
  555. ik[12] = KernelManager->AddKernel<TE, 12, INAIR, INGROUND>( );
  556. ik[10] = KernelManager->AddKernel<TE, 10, INAIR, INGROUND>( );
  557. ik[11] = KernelManager->AddKernel<TE, 11, INAIR, INGROUND>( );
  558. } else if (lays > 0 && layr == 0) {
  559. ik[12] = KernelManager->AddKernel<TE, 12, INGROUND, INAIR>( );
  560. ik[10] = KernelManager->AddKernel<TE, 10, INGROUND, INAIR>( );
  561. ik[11] = KernelManager->AddKernel<TE, 11, INGROUND, INAIR>( );
  562. } else {
  563. ik[12] = KernelManager->AddKernel<TE, 12, INGROUND, INGROUND>( );
  564. ik[10] = KernelManager->AddKernel<TE, 10, INGROUND, INGROUND>( );
  565. ik[11] = KernelManager->AddKernel<TE, 11, INGROUND, INGROUND>( );
  566. }
  567. }
  568. }
  569. if (std::abs(Pol[1]) > 0 || std::abs(Pol[0]) > 0) { // x or y grounded HED dipole
  570. switch (FieldsToCalculate) {
  571. case E:
  572. if ( lays == 0 && layr == 0) {
  573. ik[5] = KernelManager->AddKernel<TE, 5, INAIR, INAIR>( );
  574. ik[6] = KernelManager->AddKernel<TE, 6, INAIR, INAIR>( );
  575. ik[7] = KernelManager->AddKernel<TM, 7, INAIR, INAIR>( );
  576. ik[8] = KernelManager->AddKernel<TM, 8, INAIR, INAIR>( );
  577. ik[9] = KernelManager->AddKernel<TM, 9, INAIR, INAIR>( );
  578. } else if (lays == 0 && layr > 0) {
  579. ik[5] = KernelManager->AddKernel<TE, 5, INAIR, INGROUND>( );
  580. ik[6] = KernelManager->AddKernel<TE, 6, INAIR, INGROUND>( );
  581. ik[7] = KernelManager->AddKernel<TM, 7, INAIR, INGROUND>( );
  582. ik[8] = KernelManager->AddKernel<TM, 8, INAIR, INGROUND>( );
  583. ik[9] = KernelManager->AddKernel<TM, 9, INAIR, INGROUND>( );
  584. } else if (lays > 0 && layr == 0) {
  585. ik[5] = KernelManager->AddKernel<TE, 5, INGROUND, INAIR>( );
  586. ik[6] = KernelManager->AddKernel<TE, 6, INGROUND, INAIR>( );
  587. ik[7] = KernelManager->AddKernel<TM, 7, INGROUND, INAIR>( );
  588. ik[8] = KernelManager->AddKernel<TM, 8, INGROUND, INAIR>( );
  589. ik[9] = KernelManager->AddKernel<TM, 9, INGROUND, INAIR>( );
  590. } else {
  591. ik[5] = KernelManager->AddKernel<TE, 5, INGROUND, INGROUND>( );
  592. ik[6] = KernelManager->AddKernel<TE, 6, INGROUND, INGROUND>( );
  593. ik[7] = KernelManager->AddKernel<TM, 7, INGROUND, INGROUND>( );
  594. ik[8] = KernelManager->AddKernel<TM, 8, INGROUND, INGROUND>( );
  595. ik[9] = KernelManager->AddKernel<TM, 9, INGROUND, INGROUND>( );
  596. }
  597. break;
  598. case H:
  599. if ( lays == 0 && layr == 0) {
  600. ik[0] = KernelManager->AddKernel<TE, 0, INAIR, INAIR>( );
  601. ik[1] = KernelManager->AddKernel<TE, 1, INAIR, INAIR>( );
  602. ik[4] = KernelManager->AddKernel<TE, 4, INAIR, INAIR>( );
  603. ik[2] = KernelManager->AddKernel<TM, 2, INAIR, INAIR>( );
  604. ik[3] = KernelManager->AddKernel<TM, 3, INAIR, INAIR>( );
  605. } else if (lays == 0 && layr > 0) {
  606. ik[0] = KernelManager->AddKernel<TE, 0, INAIR, INGROUND>( );
  607. ik[1] = KernelManager->AddKernel<TE, 1, INAIR, INGROUND>( );
  608. ik[4] = KernelManager->AddKernel<TE, 4, INAIR, INGROUND>( );
  609. ik[2] = KernelManager->AddKernel<TM, 2, INAIR, INGROUND>( );
  610. ik[3] = KernelManager->AddKernel<TM, 3, INAIR, INGROUND>( );
  611. } else if (lays > 0 && layr == 0) {
  612. ik[0] = KernelManager->AddKernel<TE, 0, INGROUND, INAIR>( );
  613. ik[1] = KernelManager->AddKernel<TE, 1, INGROUND, INAIR>( );
  614. ik[4] = KernelManager->AddKernel<TE, 4, INGROUND, INAIR>( );
  615. ik[2] = KernelManager->AddKernel<TM, 2, INGROUND, INAIR>( );
  616. ik[3] = KernelManager->AddKernel<TM, 3, INGROUND, INAIR>( );
  617. } else {
  618. ik[0] = KernelManager->AddKernel<TE, 0, INGROUND, INGROUND>( );
  619. ik[1] = KernelManager->AddKernel<TE, 1, INGROUND, INGROUND>( );
  620. ik[4] = KernelManager->AddKernel<TE, 4, INGROUND, INGROUND>( );
  621. ik[2] = KernelManager->AddKernel<TM, 2, INGROUND, INGROUND>( );
  622. ik[3] = KernelManager->AddKernel<TM, 3, INGROUND, INGROUND>( );
  623. }
  624. break;
  625. case BOTH:
  626. if ( lays == 0 && layr == 0) {
  627. ik[5] = KernelManager->AddKernel<TE, 5, INAIR, INAIR>( );
  628. ik[6] = KernelManager->AddKernel<TE, 6, INAIR, INAIR>( );
  629. ik[7] = KernelManager->AddKernel<TM, 7, INAIR, INAIR>( );
  630. ik[8] = KernelManager->AddKernel<TM, 8, INAIR, INAIR>( );
  631. ik[9] = KernelManager->AddKernel<TM, 9, INAIR, INAIR>( );
  632. ik[0] = KernelManager->AddKernel<TE, 0, INAIR, INAIR>( );
  633. ik[1] = KernelManager->AddKernel<TE, 1, INAIR, INAIR>( );
  634. ik[4] = KernelManager->AddKernel<TE, 4, INAIR, INAIR>( );
  635. ik[2] = KernelManager->AddKernel<TM, 2, INAIR, INAIR>( );
  636. ik[3] = KernelManager->AddKernel<TM, 3, INAIR, INAIR>( );
  637. } else if (lays == 0 && layr > 0) {
  638. ik[5] = KernelManager->AddKernel<TE, 5, INAIR, INGROUND>( );
  639. ik[6] = KernelManager->AddKernel<TE, 6, INAIR, INGROUND>( );
  640. ik[7] = KernelManager->AddKernel<TM, 7, INAIR, INGROUND>( );
  641. ik[8] = KernelManager->AddKernel<TM, 8, INAIR, INGROUND>( );
  642. ik[9] = KernelManager->AddKernel<TM, 9, INAIR, INGROUND>( );
  643. ik[0] = KernelManager->AddKernel<TE, 0, INAIR, INGROUND>( );
  644. ik[1] = KernelManager->AddKernel<TE, 1, INAIR, INGROUND>( );
  645. ik[4] = KernelManager->AddKernel<TE, 4, INAIR, INGROUND>( );
  646. ik[2] = KernelManager->AddKernel<TM, 2, INAIR, INGROUND>( );
  647. ik[3] = KernelManager->AddKernel<TM, 3, INAIR, INGROUND>( );
  648. } else {
  649. ik[5] = KernelManager->AddKernel<TE, 5, INGROUND, INGROUND>( );
  650. ik[6] = KernelManager->AddKernel<TE, 6, INGROUND, INGROUND>( );
  651. ik[7] = KernelManager->AddKernel<TM, 7, INGROUND, INGROUND>( );
  652. ik[8] = KernelManager->AddKernel<TM, 8, INGROUND, INGROUND>( );
  653. ik[9] = KernelManager->AddKernel<TM, 9, INGROUND, INGROUND>( );
  654. ik[0] = KernelManager->AddKernel<TE, 0, INGROUND, INGROUND>( );
  655. ik[1] = KernelManager->AddKernel<TE, 1, INGROUND, INGROUND>( );
  656. ik[4] = KernelManager->AddKernel<TE, 4, INGROUND, INGROUND>( );
  657. ik[2] = KernelManager->AddKernel<TM, 2, INGROUND, INGROUND>( );
  658. ik[3] = KernelManager->AddKernel<TM, 3, INGROUND, INGROUND>( );
  659. }
  660. break;
  661. }
  662. }
  663. break;
  664. default:
  665. std::cerr << "Dipole type incorrect, in dipolesource.cpp";
  666. exit(EXIT_FAILURE);
  667. }
  668. }
  669. void DipoleSource::UpdateFields( const int& ifreq, HankelTransform* Hankel, const Real& wavef) {
  670. Vector3r Pol = Phat;
  671. switch (Type) {
  672. case (GROUNDEDELECTRICDIPOLE):
  673. //Hankel->ComputeRelated(rho, KernelManager);
  674. if (std::abs(Pol[2]) > 0) { // z dipole
  675. switch(FieldsToCalculate) {
  676. case E:
  677. f(10) = Hankel->Zgauss(10, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[10])) / KernelManager->GetRAWKernel(ik[10])->GetYm();
  678. f(11) = Hankel->Zgauss(11, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[11])) / KernelManager->GetRAWKernel(ik[11])->GetYm();
  679. this->Receivers->AppendEfield(ifreq, irec,
  680. -Pol[2]*QPI*cp*f(10)*Moment,
  681. -Pol[2]*QPI*sp*f(10)*Moment,
  682. Pol[2]*QPI*f(11)*Moment);
  683. break;
  684. case H:
  685. f(12) = Hankel->Zgauss(12, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[12]));
  686. this->Receivers->AppendHfield(ifreq, irec,
  687. -Pol[2]*QPI*sp*f(12)*Moment,
  688. Pol[2]*QPI*cp*f(12)*Moment,
  689. 0. );
  690. break;
  691. case BOTH:
  692. f(10) = Hankel->Zgauss(10, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[10])) / KernelManager->GetRAWKernel(ik[10])->GetYm();
  693. f(11) = Hankel->Zgauss(11, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[11])) / KernelManager->GetRAWKernel(ik[11])->GetYm();
  694. this->Receivers->AppendEfield(ifreq, irec,
  695. -Pol[2]*QPI*cp*f(10)*Moment,
  696. -Pol[2]*QPI*sp*f(10)*Moment,
  697. Pol[2]*QPI*f(11)*Moment );
  698. f(12) = Hankel->Zgauss(12, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[12]));
  699. this->Receivers->AppendHfield(ifreq, irec,
  700. -Pol[2]*QPI*sp*f(12)*Moment,
  701. Pol[2]*QPI*cp*f(12)*Moment,
  702. 0. );
  703. } // Fields to calculate Z polarity Electric dipole
  704. }
  705. if (std::abs(Pol[1]) > 0 || std::abs(Pol[0]) > 0) { // x or y dipole
  706. switch(FieldsToCalculate) {
  707. case E:
  708. f(2) = Hankel->Zgauss(2, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[2])) * KernelManager->GetRAWKernel(ik[2])->GetZs();
  709. f(3) = Hankel->Zgauss(3, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[3])) * KernelManager->GetRAWKernel(ik[3])->GetZs();
  710. f(0) = Hankel->Zgauss(0, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[0])) / KernelManager->GetRAWKernel(ik[0])->GetYm();
  711. f(1) = Hankel->Zgauss(1, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[1])) / KernelManager->GetRAWKernel(ik[1])->GetYm();
  712. f(4) = Hankel->Zgauss(4, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[4])) / KernelManager->GetRAWKernel(ik[4])->GetYm();
  713. if (std::abs(Pol[1]) > 0) {
  714. this->Receivers->AppendEfield(ifreq, irec,
  715. Pol[1]*Moment*QPI*scp*((f(0)-(Real)(2.)*f(1)/rho)+(f(2)-(Real)(2.)*f(3)/rho)),
  716. Pol[1]*Moment*QPI*((sps*f(0)+c2p*f(1)/rho)-(cps*f(2)-c2p*f(3)/rho)),
  717. Pol[1]*Moment*QPI*sp*f(4));
  718. }
  719. if (std::abs(Pol[0]) > 0) {
  720. this->Receivers->AppendEfield(ifreq, irec,
  721. Pol[0]*Moment*QPI*((cps*f(0)-c2p*f(1)/rho)-(sps*f(2)+c2p*f(3)/rho)),
  722. Pol[0]*Moment*QPI*scp*((f(0)-(Real)(2.)*f(1)/rho)+(f(2)-(Real)(2.)*f(3)/rho)),
  723. Pol[0]*Moment*QPI*cp*f(4) );
  724. }
  725. break;
  726. case H:
  727. f(5) = Hankel->Zgauss(5, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[5]));
  728. f(6) = Hankel->Zgauss(6, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[6]));
  729. f(7) = Hankel->Zgauss(7, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[7]))*KernelManager->GetRAWKernel(ik[7])->GetZs()/KernelManager->GetRAWKernel(ik[7])->GetZm();
  730. f(8) = Hankel->Zgauss(8, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[8]))*KernelManager->GetRAWKernel(ik[8])->GetZs()/KernelManager->GetRAWKernel(ik[8])->GetZm();
  731. f(9) = Hankel->Zgauss(9, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[9]))*KernelManager->GetRAWKernel(ik[9])->GetZs()/KernelManager->GetRAWKernel(ik[9])->GetZm();
  732. if (std::abs(Pol[1]) > 0) {
  733. this->Receivers->AppendHfield(ifreq, irec,
  734. Pol[1]*QPI*(sps*f(5)+c2p*f(6)/rho-cps*f(7)+c2p*f(8)/rho)*Moment,
  735. Pol[1]*QPI*scp*(-f(5)+(Real)(2.)*f(6)/rho-f(7)+(Real)(2.)*f(8)/rho)*Moment,
  736. -Pol[1]*QPI*cp*f(9)*Moment );
  737. }
  738. if (std::abs(Pol[0]) > 0) {
  739. this->Receivers->AppendHfield(ifreq, irec,
  740. Pol[0]*Moment*QPI*scp*(f(5)-(Real)(2.)*f(6)/rho+f(7)-(Real)(2.)*f(8)/rho),
  741. Pol[0]*Moment*QPI*(-cps*f(5)+c2p*f(6)/rho+sps*f(7)+c2p*f(8)/rho),
  742. Pol[0]*Moment*QPI*sp*f(9) );
  743. }
  744. break;
  745. case BOTH:
  746. f(0) = Hankel->Zgauss(0, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[0])) / KernelManager->GetRAWKernel(ik[0])->GetYm();
  747. f(1) = Hankel->Zgauss(1, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[1])) / KernelManager->GetRAWKernel(ik[1])->GetYm();
  748. f(4) = Hankel->Zgauss(4, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[4])) / KernelManager->GetRAWKernel(ik[4])->GetYm();
  749. f(2) = Hankel->Zgauss(2, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[2])) * KernelManager->GetRAWKernel(ik[2])->GetZs();
  750. f(3) = Hankel->Zgauss(3, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[3])) * KernelManager->GetRAWKernel(ik[3])->GetZs();
  751. f(5) = Hankel->Zgauss(5, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[5]));
  752. f(6) = Hankel->Zgauss(6, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[6]));
  753. f(7) = Hankel->Zgauss(7, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[7]))*KernelManager->GetRAWKernel(ik[7])->GetZs()/KernelManager->GetRAWKernel(ik[7])->GetZm();
  754. f(8) = Hankel->Zgauss(8, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[8]))*KernelManager->GetRAWKernel(ik[8])->GetZs()/KernelManager->GetRAWKernel(ik[8])->GetZm();
  755. f(9) = Hankel->Zgauss(9, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[9]))*KernelManager->GetRAWKernel(ik[9])->GetZs()/KernelManager->GetRAWKernel(ik[9])->GetZm();
  756. if (std::abs(Pol[1]) > 0) {
  757. this->Receivers->AppendEfield(ifreq, irec,
  758. Pol[1]*QPI*scp*((f(0)-(Real)(2.)*f(1)/rho)+(f(2)-(Real)(2.)*f(3)/rho))*Moment ,
  759. Pol[1]*QPI*((sps*f(0)+c2p*f(1)/rho)-(cps*f(2)-c2p*f(3)/rho))*Moment,
  760. Pol[1]*QPI*sp*f(4)*Moment);
  761. this->Receivers->AppendHfield(ifreq, irec,
  762. Pol[1]*QPI*(sps*f(5)+c2p*f(6)/rho-cps*f(7)+c2p*f(8)/rho)*Moment,
  763. Pol[1]*QPI*scp*(-f(5)+(Real)(2.)*f(6)/rho-f(7)+(Real)(2.)*f(8)/rho)*Moment,
  764. -Pol[1]*QPI*cp*f(9)*Moment );
  765. }
  766. if (std::abs(Pol[0]) > 0) {
  767. this->Receivers->AppendEfield(ifreq, irec,
  768. Pol[0]*Moment*QPI*((cps*f(0)-c2p*f(1)/rho)-(sps*f(2)+c2p*f(3)/rho)),
  769. Pol[0]*Moment*QPI*scp*((f(0)-(Real)(2.)*f(1)/rho)+(f(2)-(Real)(2.)*f(3)/rho)),
  770. Pol[0]*Moment*QPI*cp*f(4) );
  771. this->Receivers->AppendHfield(ifreq, irec,
  772. Pol[0]*Moment*QPI*scp*(f(5)-(Real)(2.)*f(6)/rho+f(7)-(Real)(2.)*f(8)/rho),
  773. Pol[0]*Moment*QPI*(-cps*f(5)+c2p*f(6)/rho+sps*f(7)+c2p*f(8)/rho),
  774. Pol[0]*Moment*QPI*sp*f(9) );
  775. }
  776. break;
  777. }
  778. }
  779. break; // GROUNDEDELECTRICDIPOLE
  780. case UNGROUNDEDELECTRICDIPOLE:
  781. if (std::abs(Pol[2]) > 0) { // z dipole
  782. switch(FieldsToCalculate) {
  783. case E:
  784. f(10) = 0;
  785. f(11) = Hankel->Zgauss(11, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[11])) / KernelManager->GetRAWKernel(ik[11])->GetYm();
  786. this->Receivers->AppendEfield(ifreq, irec,
  787. -Pol[2]*QPI*cp*f(10)*Moment,
  788. -Pol[2]*QPI*sp*f(10)*Moment,
  789. Pol[2]*QPI*f(11)*Moment);
  790. break;
  791. case H:
  792. f(12) = Hankel->Zgauss(12, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[12]));
  793. this->Receivers->AppendHfield(ifreq, irec,
  794. -Pol[2]*QPI*sp*f(12)*Moment,
  795. Pol[2]*QPI*cp*f(12)*Moment,
  796. 0. );
  797. break;
  798. case BOTH:
  799. f(10) = 0;
  800. f(11) = Hankel->Zgauss(11, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[11])) / KernelManager->GetRAWKernel(ik[11])->GetYm();
  801. this->Receivers->AppendEfield(ifreq, irec,
  802. -Pol[2]*QPI*cp*f(10)*Moment,
  803. -Pol[2]*QPI*sp*f(10)*Moment,
  804. Pol[2]*QPI*f(11)*Moment );
  805. f(12) = Hankel->Zgauss(12, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[12]));
  806. this->Receivers->AppendHfield(ifreq, irec,
  807. -Pol[2]*QPI*sp*f(12)*Moment,
  808. Pol[2]*QPI*cp*f(12)*Moment,
  809. 0. );
  810. } // Fields to calculate Z polarity Electric dipole
  811. }
  812. if (std::abs(Pol[1]) > 0 || std::abs(Pol[0]) > 0) { // x or y dipole
  813. switch(FieldsToCalculate) {
  814. case E:
  815. f(0) = 0;
  816. f(1) = 0;
  817. f(2) = Hankel->Zgauss(2, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[2])) * KernelManager->GetRAWKernel(ik[2])->GetZs();
  818. f(3) = Hankel->Zgauss(3, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[3])) * KernelManager->GetRAWKernel(ik[3])->GetZs();
  819. f(4) = 0;
  820. if (std::abs(Pol[1]) > 0) {
  821. this->Receivers->AppendEfield(ifreq, irec,
  822. Pol[1]*QPI*scp*((f(0)-(Real)(2.)*f(1)/rho)+(f(2)-(Real)(2.)*f(3)/rho))*Moment,
  823. Pol[1]*QPI*((sps*f(0)+c2p*f(1)/rho)-(cps*f(2)-c2p*f(3)/rho))*Moment,
  824. Pol[1]*QPI*sp*f(4)*Moment);
  825. }
  826. if (std::abs(Pol[0]) > 0) {
  827. this->Receivers->AppendEfield(ifreq, irec,
  828. Pol[0]*Moment*QPI*((cps*f(0)-c2p*f(1)/rho)-(sps*f(2)+c2p*f(3)/rho)),
  829. Pol[0]*Moment*QPI*scp*((f(0)-(Real)(2.)*f(1)/rho)+(f(2)-(Real)(2.)*f(3)/rho)),
  830. Pol[0]*Moment*QPI*cp*f(4) );
  831. }
  832. break;
  833. case H:
  834. f(5) = Hankel->Zgauss(5, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[5]));
  835. f(6) = Hankel->Zgauss(6, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[6]));
  836. f(7) = Hankel->Zgauss(7, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[7]))*KernelManager->GetRAWKernel(ik[7])->GetZs()/KernelManager->GetRAWKernel(ik[7])->GetZm();
  837. f(8) = Hankel->Zgauss(8, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[8]))*KernelManager->GetRAWKernel(ik[8])->GetZs()/KernelManager->GetRAWKernel(ik[8])->GetZm();
  838. f(9) = Hankel->Zgauss(9, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[9]))*KernelManager->GetRAWKernel(ik[9])->GetZs()/KernelManager->GetRAWKernel(ik[9])->GetZm();
  839. if (std::abs(Pol[1]) > 0) {
  840. this->Receivers->AppendHfield(ifreq, irec,
  841. Pol[1]*QPI*(sps*f(5)+c2p*f(6)/rho-cps*f(7)+c2p*f(8)/rho)*Moment,
  842. Pol[1]*QPI*scp*(-f(5)+(Real)(2.)*f(6)/rho-f(7)+(Real)(2.)*f(8)/rho)*Moment,
  843. -Pol[1]*QPI*cp*f(9)*Moment );
  844. // Analytic whole space solution could go here
  845. }
  846. if (std::abs(Pol[0]) > 0) {
  847. this->Receivers->AppendHfield(ifreq, irec,
  848. Pol[0]*Moment*QPI*scp*(f(5)-(Real)(2.)*f(6)/rho+f(7)-(Real)(2.)*f(8)/rho),
  849. Pol[0]*Moment*QPI*(-cps*f(5)+c2p*f(6)/rho+sps*f(7)+c2p*f(8)/rho),
  850. Pol[0]*Moment*QPI*sp*f(9) );
  851. // Analytic whole space solution
  852. }
  853. break;
  854. case BOTH:
  855. f(0) = 0;
  856. f(1) = 0;
  857. f(4) = 0;
  858. f(2) = Hankel->Zgauss(2, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[2])) * KernelManager->GetRAWKernel(0)->GetZs();
  859. f(3) = Hankel->Zgauss(3, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[3])) * KernelManager->GetRAWKernel(1)->GetZs();
  860. f(5) = Hankel->Zgauss(5, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[5]));
  861. f(6) = Hankel->Zgauss(6, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[6]));
  862. f(7) = Hankel->Zgauss(7, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[7]))*KernelManager->GetRAWKernel(ik[7])->GetZs()/KernelManager->GetRAWKernel(ik[7])->GetZm();
  863. f(8) = Hankel->Zgauss(8, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[8]))*KernelManager->GetRAWKernel(ik[8])->GetZs()/KernelManager->GetRAWKernel(ik[8])->GetZm();
  864. f(9) = Hankel->Zgauss(9, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[9]))*KernelManager->GetRAWKernel(ik[9])->GetZs()/KernelManager->GetRAWKernel(ik[9])->GetZm();
  865. if (std::abs(Pol[1]) > 0) {
  866. this->Receivers->AppendEfield(ifreq, irec,
  867. Pol[1]*QPI*scp*((f(0)-(Real)(2.)*f(1)/rho)+(f(2)-(Real)(2.)*f(3)/rho))*Moment ,
  868. Pol[1]*QPI*((sps*f(0)+c2p*f(1)/rho)-(cps*f(2)-c2p*f(3)/rho))*Moment,
  869. Pol[1]*QPI*sp*f(4)*Moment);
  870. this->Receivers->AppendHfield(ifreq, irec,
  871. Pol[1]*QPI*(sps*f(5)+c2p*f(6)/rho-cps*f(7)+c2p*f(8)/rho)*Moment,
  872. Pol[1]*QPI*scp*(-f(5)+(Real)(2.)*f(6)/rho-f(7)+(Real)(2.)*f(8)/rho)*Moment,
  873. -Pol[1]*QPI*cp*f(9)*Moment );
  874. }
  875. if (std::abs(Pol[0]) > 0) {
  876. this->Receivers->AppendEfield(ifreq, irec,
  877. Pol[0]*Moment*QPI*((cps*f(0)-c2p*f(1)/rho)-(sps*f(2)+c2p*f(3)/rho)),
  878. Pol[0]*Moment*QPI*scp*((f(0)-(Real)(2.)*f(1)/rho)+(f(2)-(Real)(2.)*f(3)/rho)),
  879. Pol[0]*Moment*QPI*cp*f(4) );
  880. this->Receivers->AppendHfield(ifreq, irec,
  881. Pol[0]*Moment*QPI*scp*(f(5)-(Real)(2.)*f(6)/rho+f(7)-(Real)(2.)*f(8)/rho),
  882. Pol[0]*Moment*QPI*(-cps*f(5)+c2p*f(6)/rho+sps*f(7)+c2p*f(8)/rho),
  883. Pol[0]*Moment*QPI*sp*f(9) );
  884. }
  885. break;
  886. }
  887. }
  888. break; // UNGROUNDEDELECTRICDIPOLE
  889. case MAGNETICDIPOLE:
  890. //Hankel->ComputeRelated(rho, KernelManager);
  891. if (std::abs(Pol[2]) > 0) { // z dipole
  892. switch(FieldsToCalculate) {
  893. case E:
  894. f(12)=Hankel->Zgauss(12, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[12]))*KernelManager->GetRAWKernel(ik[12])->GetZs();
  895. this->Receivers->AppendEfield(ifreq, irec,
  896. Pol[2]*Moment*QPI*sp*f(12),
  897. -Pol[2]*Moment*QPI*cp*f(12),
  898. 0);
  899. break;
  900. case H:
  901. f(10)=Hankel->Zgauss(10, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[10]))*KernelManager->GetRAWKernel(ik[10])->GetZs()/KernelManager->GetRAWKernel(ik[10])->GetZm();
  902. f(11)=Hankel->Zgauss(11, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[11]))*KernelManager->GetRAWKernel(ik[11])->GetZs()/KernelManager->GetRAWKernel(ik[11])->GetZm();
  903. this->Receivers->AppendHfield(ifreq, irec,
  904. -Pol[2]*Moment*QPI*cp*f(10),
  905. -Pol[2]*Moment*QPI*sp*f(10),
  906. Pol[2]*Moment*QPI*f(11) );
  907. break;
  908. case BOTH:
  909. f(12)=Hankel->Zgauss(12, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[12]))*KernelManager->GetRAWKernel(ik[12])->GetZs();
  910. f(10)=Hankel->Zgauss(10, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[10]))*KernelManager->GetRAWKernel(ik[10])->GetZs()/KernelManager->GetRAWKernel(ik[10])->GetZm();
  911. f(11)=Hankel->Zgauss(11, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[11]))*KernelManager->GetRAWKernel(ik[11])->GetZs()/KernelManager->GetRAWKernel(ik[11])->GetZm();
  912. this->Receivers->AppendEfield(ifreq, irec,
  913. Pol[2]*Moment*QPI*sp*f(12),
  914. -Pol[2]*Moment*QPI*cp*f(12),
  915. 0);
  916. this->Receivers->AppendHfield(ifreq, irec,
  917. -Pol[2]*Moment*QPI*cp*f(10),
  918. -Pol[2]*Moment*QPI*sp*f(10),
  919. Pol[2]*Moment*QPI*f(11) );
  920. break;
  921. }
  922. }
  923. if (std::abs(Pol[1]) > 0 || std::abs(Pol[0]) > 0) { // x or y grounded HED dipole
  924. switch (FieldsToCalculate) {
  925. case E:
  926. f(5) = Hankel->Zgauss(5, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[5]))*KernelManager->GetRAWKernel(ik[5])->GetZs();
  927. f(6) = Hankel->Zgauss(6, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[6]))*KernelManager->GetRAWKernel(ik[6])->GetZs();
  928. f(7) = Hankel->Zgauss(7, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[7]))*KernelManager->GetRAWKernel(ik[7])->GetKs()/KernelManager->GetRAWKernel(ik[7])->GetYm();
  929. f(8) = Hankel->Zgauss(8, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[8]))*KernelManager->GetRAWKernel(ik[8])->GetKs()/KernelManager->GetRAWKernel(ik[8])->GetYm();
  930. f(9) = Hankel->Zgauss(9, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[9]))*KernelManager->GetRAWKernel(ik[9])->GetKs()/KernelManager->GetRAWKernel(ik[9])->GetYm();
  931. if (std::abs(Pol[0]) > 0) {
  932. this->Receivers->AppendEfield(ifreq, irec,
  933. Pol[0]*Moment*QPI*scp*((-f(5)+(Real)(2.)*f(6)/rho)+(f(7)-(Real)(2.)*f(8)/rho)),
  934. Pol[0]*Moment*QPI*((cps*f(5)-c2p*f(6)/rho)+(sps*f(7)+c2p*f(8)/rho)),
  935. Pol[0]*Moment*QPI*sp*f(9));
  936. }
  937. if (std::abs(Pol[1]) > 0) {
  938. this->Receivers->AppendEfield(ifreq, irec,
  939. Pol[1]*Moment*QPI*(-(sps*f(5)+c2p*f(6)/rho)-(cps*f(7)-c2p*f(8)/rho)),
  940. Pol[1]*Moment*QPI*scp*((f(5)-(Real)(2.)*f(6)/rho)-(f(7)-(Real)(2.)*f(8)/rho)),
  941. -Pol[1]*Moment*QPI*cp*f(9) );
  942. }
  943. break;
  944. case H:
  945. f(0) = Hankel->Zgauss(0, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[0]))*KernelManager->GetRAWKernel(ik[0])->GetZs()/KernelManager->GetRAWKernel(ik[0])->GetZm();
  946. f(1) = Hankel->Zgauss(1, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[1]))*KernelManager->GetRAWKernel(ik[1])->GetZs()/KernelManager->GetRAWKernel(ik[1])->GetZm();
  947. f(4) = Hankel->Zgauss(4, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[4]))*KernelManager->GetRAWKernel(ik[4])->GetZs()/KernelManager->GetRAWKernel(ik[4])->GetZm();
  948. f(2) = Hankel->Zgauss(2, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[2]))*KernelManager->GetRAWKernel(ik[2])->GetKs();
  949. f(3) = Hankel->Zgauss(3, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[3]))*KernelManager->GetRAWKernel(ik[3])->GetKs();
  950. if (std::abs(Pol[0]) > 0) {
  951. this->Receivers->AppendHfield(ifreq, irec,
  952. Pol[0]*Moment*QPI*(cps*f(0)-c2p*f(1)/rho+(sps*f(2)+c2p*f(3)/rho)),
  953. Pol[0]*Moment*QPI*scp*(f(0)-(Real)(2.)*f(1)/rho-(f(2)-(Real)(2.)*f(3)/rho)),
  954. Pol[0]*Moment*QPI*cp*f(4) );
  955. }
  956. if (std::abs(Pol[1]) > 0) {
  957. this->Receivers->AppendHfield(ifreq, irec,
  958. Pol[1]*Moment*QPI*scp*(f(0)-(Real)(2.)*f(1)/rho-(f(2)-(Real)(2.)*f(3)/rho)),
  959. Pol[1]*Moment*QPI*(sps*f(0)+c2p*f(1)/rho+(cps*f(2)-c2p*f(3)/rho)),
  960. Pol[1]*Moment*QPI*sp*f(4));
  961. }
  962. break;
  963. case BOTH:
  964. f(5) = Hankel->Zgauss(5, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[5]))*KernelManager->GetRAWKernel(ik[5])->GetZs();
  965. f(6) = Hankel->Zgauss(6, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[6]))*KernelManager->GetRAWKernel(ik[6])->GetZs();
  966. f(7) = Hankel->Zgauss(7, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[7]))*KernelManager->GetRAWKernel(ik[7])->GetKs()/KernelManager->GetRAWKernel(ik[7])->GetYm();
  967. f(8) = Hankel->Zgauss(8, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[8]))*KernelManager->GetRAWKernel(ik[8])->GetKs()/KernelManager->GetRAWKernel(ik[8])->GetYm();
  968. f(9) = Hankel->Zgauss(9, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[9]))*KernelManager->GetRAWKernel(ik[9])->GetKs()/KernelManager->GetRAWKernel(ik[9])->GetYm();
  969. f(0) = Hankel->Zgauss(0, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[0]))*KernelManager->GetRAWKernel(ik[0])->GetZs()/KernelManager->GetRAWKernel(ik[0])->GetZm();
  970. f(1) = Hankel->Zgauss(1, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[1]))*KernelManager->GetRAWKernel(ik[1])->GetZs()/KernelManager->GetRAWKernel(ik[1])->GetZm();
  971. f(4) = Hankel->Zgauss(4, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[4]))*KernelManager->GetRAWKernel(ik[4])->GetZs()/KernelManager->GetRAWKernel(ik[4])->GetZm();
  972. f(2) = Hankel->Zgauss(2, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[2]))*KernelManager->GetRAWKernel(ik[2])->GetKs();
  973. f(3) = Hankel->Zgauss(3, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[3]))*KernelManager->GetRAWKernel(ik[3])->GetKs();
  974. if (std::abs(Pol[0]) > 0) {
  975. this->Receivers->AppendEfield(ifreq, irec,
  976. Pol[0]*Moment*QPI*scp*((-f(5)+(Real)(2.)*f(6)/rho)+(f(7)-(Real)(2.)*f(8)/rho)),
  977. Pol[0]*Moment*QPI*((cps*f(5)-c2p*f(6)/rho)+(sps*f(7)+c2p*f(8)/rho)),
  978. Pol[0]*Moment*QPI*sp*f(9));
  979. this->Receivers->AppendHfield(ifreq, irec,
  980. Pol[0]*Moment*QPI*(cps*f(0)-c2p*f(1)/rho+(sps*f(2)+c2p*f(3)/rho)),
  981. Pol[0]*Moment*QPI*scp*(f(0)-(Real)(2.)*f(1)/rho-(f(2)-(Real)(2.)*f(3)/rho)),
  982. Pol[0]*Moment*QPI*cp*f(4) );
  983. }
  984. if (std::abs(Pol[1]) > 0) {
  985. this->Receivers->AppendEfield(ifreq, irec,
  986. Pol[1]*Moment*QPI*(-(sps*f(5)+c2p*f(6)/rho)-(cps*f(7)-c2p*f(8)/rho)),
  987. Pol[1]*Moment*QPI*scp*((f(5)-(Real)(2.)*f(6)/rho)-(f(7)-(Real)(2.)*f(8)/rho)),
  988. -Pol[1]*Moment*QPI*cp*f(9) );
  989. this->Receivers->AppendHfield(ifreq, irec,
  990. Pol[1]*Moment*QPI*scp*(f(0)-(Real)(2.)*f(1)/rho-(f(2)-(Real)(2.)*f(3)/rho)),
  991. Pol[1]*Moment*QPI*(sps*f(0)+c2p*f(1)/rho+(cps*f(2)-c2p*f(3)/rho)),
  992. Pol[1]*Moment*QPI*sp*f(4));
  993. }
  994. break;
  995. }
  996. }
  997. break;
  998. case NOSOURCETYPE:
  999. throw NonValidDipoleType(this);
  1000. } // Source Type Switch
  1001. }
  1002. // ==================== INQUIRY ======================
  1003. std::shared_ptr<KernelEM1DManager> DipoleSource::GetKernelManager() {
  1004. return KernelManager;
  1005. }
  1006. Vector3r DipoleSource::GetLocation() {
  1007. return this->Location;
  1008. }
  1009. #ifdef LEMMAUSEVTK
  1010. vtkActor* DipoleSource::GetVtkActor() {
  1011. vtkActor* vActor;
  1012. vtkLineSource* vLineSource;
  1013. vtkTubeFilter* vTube;
  1014. vtkPolyDataMapper* vMapper;
  1015. vtkRegularPolygonSource* vCircleSource;
  1016. vLineSource = vtkLineSource::New();
  1017. vTube = vtkTubeFilter::New();
  1018. vMapper = vtkPolyDataMapper::New();
  1019. vCircleSource = vtkRegularPolygonSource::New();
  1020. VectorXr M0 = Location - .5*Moment*Phat;
  1021. VectorXr M1 = Location + .5*Moment*Phat;
  1022. vActor = vtkActor::New();
  1023. switch (Type) {
  1024. case GROUNDEDELECTRICDIPOLE:
  1025. vLineSource->SetPoint1( M0(0), M0(1), M0(2));
  1026. vLineSource->SetPoint2( M1(0), M1(1), M1(2));
  1027. vTube->SetInputConnection(vLineSource->GetOutputPort());
  1028. vTube->SetRadius(.1 * std::abs(Moment));
  1029. vTube->SetNumberOfSides(6);
  1030. vTube->SetCapping(1);
  1031. vMapper->SetInputConnection(vTube->GetOutputPort());
  1032. vActor->SetMapper(vMapper);
  1033. vActor->GetProperty()->SetColor(Phat[0], Phat[1], Phat[2]);
  1034. break;
  1035. case UNGROUNDEDELECTRICDIPOLE:
  1036. vLineSource->SetPoint1( M0(0), M0(1), M0(2));
  1037. vLineSource->SetPoint2( M1(0), M1(1), M1(2));
  1038. vTube->SetInputConnection(vLineSource->GetOutputPort());
  1039. vTube->SetRadius(.1 * std::abs(Moment));
  1040. vTube->SetNumberOfSides(6);
  1041. vTube->SetCapping(1);
  1042. vMapper->SetInputConnection(vTube->GetOutputPort());
  1043. vActor->SetMapper(vMapper);
  1044. //vActor->GetProperty()->SetColor(Phat[0], Phat[1], Phat[2]);
  1045. vActor->GetProperty()->SetColor(rand()/(Real)(RAND_MAX), rand()/(Real)(RAND_MAX), rand()/(Real)(RAND_MAX));
  1046. vActor->GetProperty()->SetOpacity(1.);
  1047. break;
  1048. case MAGNETICDIPOLE:
  1049. vCircleSource->SetCenter(Location(0), Location(1),
  1050. Location(2));
  1051. vCircleSource->SetNumberOfSides(360);
  1052. vCircleSource->SetNormal(Phat[0], Phat[1], Phat[2]);
  1053. vCircleSource->SetRadius(0.2); // .2 m radius
  1054. vCircleSource->SetGeneratePolygon(false);
  1055. vCircleSource->SetGeneratePolyline(true);
  1056. vCircleSource->Update();
  1057. vTube->SetInputConnection(vCircleSource->GetOutputPort());
  1058. //vTube->SetRadius( max((float)(*xCoords->GetTuple(nx)),
  1059. // (float)(*yCoords->GetTuple(ny))) / 100);
  1060. vTube->SetRadius(.1*std::abs(Moment));
  1061. vTube->SetNumberOfSides(6);
  1062. vTube->SetCapping(1);
  1063. vMapper->SetInputConnection(vTube->GetOutputPort());
  1064. vActor->SetMapper(vMapper);
  1065. vActor->GetProperty()->SetColor(.9,.2,.9);
  1066. break;
  1067. default:
  1068. throw NonValidDipoleType();
  1069. }
  1070. vLineSource->Delete();
  1071. vCircleSource->Delete();
  1072. vTube->Delete();
  1073. vMapper->Delete();
  1074. return vActor;
  1075. }
  1076. #endif
  1077. Real DipoleSource::GetLocation(const int& coordinate) {
  1078. switch (coordinate) {
  1079. case (0):
  1080. return this->Location.x();
  1081. //break; // implicit
  1082. case (1):
  1083. return this->Location.y();
  1084. //break; // implicit
  1085. case (2):
  1086. return this->Location.z();
  1087. //break; // implicit
  1088. default:
  1089. throw NonValidLocationCoordinate( );
  1090. }
  1091. }
  1092. DIPOLESOURCETYPE DipoleSource::GetDipoleSourceType() {
  1093. return this->Type;
  1094. }
  1095. //DipoleSourcePolarisation DipoleSource::GetDipoleSourcePolarisation() {
  1096. // return this->Polarisation;
  1097. //}
  1098. Real DipoleSource::GetAngularFrequency(const int& ifreq) {
  1099. return 2.*PI*this->Freqs(ifreq);
  1100. }
  1101. Real DipoleSource::GetFrequency(const int& ifreq) {
  1102. return this->Freqs(ifreq);
  1103. }
  1104. VectorXr DipoleSource::GetFrequencies( ) {
  1105. return this->Freqs;
  1106. }
  1107. Real DipoleSource::GetPhase() {
  1108. return this->Phase;
  1109. }
  1110. Real DipoleSource::GetMoment() {
  1111. return this->Moment;
  1112. }
  1113. int DipoleSource::GetNumberOfFrequencies() {
  1114. return (int)(this->Freqs.size());
  1115. }
  1116. void DipoleSource::SetNumberOfFrequencies(const int &nfreq){
  1117. Freqs.resize(nfreq);
  1118. Freqs.setZero();
  1119. }
  1120. void DipoleSource::SetFrequency(const int &ifreq, const Real &freq){
  1121. Freqs(ifreq) = freq;
  1122. }
  1123. void DipoleSource::SetFrequencies(const VectorXr &freqs){
  1124. Freqs = freqs;
  1125. }
  1126. /////////////////////////////////////////////////////////////////
  1127. /////////////////////////////////////////////////////////////////
  1128. // Error classes
  1129. NullDipoleSource::NullDipoleSource() :
  1130. runtime_error( "NULL VALUED DIPOLE SOURCE") {}
  1131. NonValidDipoleTypeAssignment::NonValidDipoleTypeAssignment( ) :
  1132. runtime_error( "NON VALID DIPOLE TYPE ASSIGNMENT") { }
  1133. NonValidDipoleType::NonValidDipoleType( LemmaObject* ptr ) :
  1134. runtime_error( "NON VALID DIPOLE TYPE") {
  1135. std::cout << "Thrown by instance of "
  1136. << ptr->GetName() << std::endl;
  1137. }
  1138. NonValidDipoleType::NonValidDipoleType( ) :
  1139. runtime_error( "NON VALID DIPOLE TYPE") { }
  1140. NonValidDipolePolarity::NonValidDipolePolarity () :
  1141. runtime_error( "NON VALID DIPOLE POLARITY") { }
  1142. NonValidDipolePolarisation::NonValidDipolePolarisation( ) :
  1143. runtime_error( "NON VALID DIPOLE TYPE") { }
  1144. NonValidDipolePolarisationAssignment::
  1145. NonValidDipolePolarisationAssignment( ) :
  1146. runtime_error( "NON VALID DIPOLE POLARISATION ASSIGNMENT") { }
  1147. NonValidLocationCoordinate::NonValidLocationCoordinate( ) :
  1148. runtime_error( "NON VALID LOCATION COORDINATE REQUESTED") { }
  1149. }