Lemma is an Electromagnetics API
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

EMEarth1D.cpp 37KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893
  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 "EMEarth1D.h"
  11. #include "FieldPoints.h"
  12. #include "WireAntenna.h"
  13. #include "PolygonalWireAntenna.h"
  14. #ifdef LEMMAUSEOMP
  15. #include "omp.h"
  16. #endif
  17. namespace Lemma {
  18. std::ostream &operator << (std::ostream &stream, const EMEarth1D &ob) {
  19. stream << ob.Serialize() << "\n";
  20. return stream;
  21. }
  22. #ifdef KIHALEE_EM1D
  23. // Wrapper function for Fortran subroutine Em1D bi kihand
  24. // Returns E or H fields (SLOW)
  25. extern "C" { void em1dcall_(int &itype, // source
  26. int &ipol, // source
  27. int &nlay, // Earth
  28. int &nfreq, // source
  29. int &nfield, // Calculator
  30. int &nres, // Receivers
  31. int &jtype, // N/A
  32. int &jgamma, // Controller
  33. double &acc, // Controller
  34. double *dep, // Earth
  35. std::complex<double> *sig, // Earth
  36. double *susl, // Earth
  37. double *sush, // Earth
  38. double *sustau, // Earth
  39. double *susalp, // Earth
  40. double *eprl, // Earth
  41. double *eprh, // Earth
  42. double *eprtau, // Earth
  43. double *epralp, // Earth
  44. double &finit, // N/A
  45. double &flimit, // N/A
  46. double &dlimit, // N/A
  47. double &lfinc, // N/A
  48. double &tx, // Source
  49. double &ty, // Source
  50. double &tz, // Source
  51. double *rxx, // Receivers
  52. double *rxy, // Receivers
  53. double *rxz, // Receivers
  54. std::complex<double> *ex, // Receivers
  55. std::complex<double> *ey, // |
  56. std::complex<double> *ez, // |
  57. std::complex<double> *hx, // |
  58. std::complex<double> *hy, // V
  59. std::complex<double> *hz ); // ___
  60. }
  61. #endif
  62. // ==================== LIFECYCLE ===================================
  63. // TODO init large arrays here.
  64. EMEarth1D::EMEarth1D( const ctor_key& key ) : LemmaObject( key ),
  65. Dipole(nullptr), Earth(nullptr), Receivers(nullptr), Antenna(nullptr),
  66. FieldsToCalculate(BOTH), HankelType(ANDERSON801), icalcinner(0), icalc(0)
  67. //#ifdef HAVE_BOOST_PROGRESS
  68. // , disp(0)
  69. //#endif
  70. {
  71. }
  72. EMEarth1D::~EMEarth1D() {
  73. }
  74. std::shared_ptr<EMEarth1D> EMEarth1D::NewSP() {
  75. return std::make_shared<EMEarth1D>(ctor_key());
  76. }
  77. YAML::Node EMEarth1D::Serialize() const {
  78. YAML::Node node = LemmaObject::Serialize();
  79. node["FieldsToCalculate"] = enum2String(FieldsToCalculate);
  80. node["HankelType"] = enum2String(HankelType);
  81. //if (Dipole != nullptr) node["Dipole"] = Dipole->Serialize();
  82. if (Earth != nullptr) node["Earth"] = Earth->Serialize();
  83. //if (Receivers != nullptr) node["Receivers"] = Receivers->Serialize(); Can be huge?
  84. if (Antenna != nullptr) node["Antenna"] = Antenna->Serialize();
  85. node.SetTag( this->GetName() );
  86. return node;
  87. }
  88. //--------------------------------------------------------------------------------------
  89. // Class: EMEarth1D
  90. // Method: GetName
  91. // Description: Class identifier
  92. //--------------------------------------------------------------------------------------
  93. inline std::string EMEarth1D::GetName ( ) const {
  94. return CName;
  95. } // ----- end of method EMEarth1D::GetName -----
  96. // ==================== ACCESS ===================================
  97. void EMEarth1D::AttachDipoleSource( std::shared_ptr<DipoleSource> dipoleptr) {
  98. Dipole = dipoleptr;
  99. }
  100. void EMEarth1D::AttachLayeredEarthEM( std::shared_ptr<LayeredEarthEM> earthptr) {
  101. Earth = earthptr;
  102. }
  103. void EMEarth1D::AttachFieldPoints( std::shared_ptr<FieldPoints> recptr) {
  104. Receivers = recptr;
  105. if (Receivers == nullptr) {
  106. std::cout << "nullptr Receivers in emearth1d.cpp " << std::endl;
  107. return;
  108. }
  109. // This has an implicid need to first set a source before receivers, users
  110. // will not expect this. Fix
  111. if (Dipole != nullptr) {
  112. switch (FieldsToCalculate) {
  113. case E:
  114. Receivers->SetNumberOfBinsE(Dipole->GetNumberOfFrequencies());
  115. break;
  116. case H:
  117. Receivers->SetNumberOfBinsH(Dipole->GetNumberOfFrequencies());
  118. break;
  119. case BOTH:
  120. Receivers->SetNumberOfBinsE(Dipole->GetNumberOfFrequencies());
  121. Receivers->SetNumberOfBinsH(Dipole->GetNumberOfFrequencies());
  122. break;
  123. }
  124. } else if (Antenna != nullptr) {
  125. switch (FieldsToCalculate) {
  126. case E:
  127. Receivers->SetNumberOfBinsE(Antenna->GetNumberOfFrequencies());
  128. break;
  129. case H:
  130. Receivers->SetNumberOfBinsH(Antenna->GetNumberOfFrequencies());
  131. break;
  132. case BOTH:
  133. Receivers->SetNumberOfBinsE(Antenna->GetNumberOfFrequencies());
  134. Receivers->SetNumberOfBinsH(Antenna->GetNumberOfFrequencies());
  135. break;
  136. }
  137. }
  138. }
  139. void EMEarth1D::AttachWireAntenna(std::shared_ptr<WireAntenna> antennae) {
  140. this->Antenna = antennae;
  141. }
  142. void EMEarth1D::SetFieldsToCalculate(const FIELDCALCULATIONS &calc) {
  143. FieldsToCalculate = calc;
  144. }
  145. void EMEarth1D::SetHankelTransformMethod( const HANKELTRANSFORMTYPE &type) {
  146. HankelType = type;
  147. }
  148. void EMEarth1D::Query() {
  149. std::cout << "EmEarth1D::Query()" << std::endl;
  150. std::cout << "Dipole " << Dipole;
  151. if (Dipole) std::cout << *Dipole << std::endl;
  152. std::cout << "Earth " << Earth;
  153. if (Earth) std::cout << *Earth << std::endl;
  154. std::cout << "Receivers " << Earth;
  155. if (Earth) std::cout << *Receivers << std::endl;
  156. std::cout << "Antenna " << Earth;
  157. if (Antenna) std::cout << *Antenna << std::endl;
  158. std::cout << "icalc " << icalc << std::endl;
  159. std::cout << "icalcinner " << icalcinner << std::endl;
  160. }
  161. // ==================== OPERATIONS ===================================
  162. void EMEarth1D::CalculateWireAntennaFields(bool progressbar) {
  163. #ifdef HAVE_BOOST_PROGRESS
  164. boost::progress_display *disp;
  165. #endif
  166. if (Earth == nullptr) {
  167. throw NullEarth();
  168. }
  169. if (Receivers == nullptr) {
  170. throw NullReceivers();
  171. }
  172. if (Antenna == nullptr) {
  173. throw NullAntenna();
  174. }
  175. if (Dipole != nullptr) {
  176. throw DipoleSourceSpecifiedForWireAntennaCalc();
  177. }
  178. Receivers->ClearFields();
  179. // Check to make sure Receivers are set up for all calculations
  180. switch(FieldsToCalculate) {
  181. case E:
  182. if (Receivers->NumberOfBinsE != Antenna->GetNumberOfFrequencies())
  183. Receivers->SetNumberOfBinsE(Antenna->GetNumberOfFrequencies());
  184. break;
  185. case H:
  186. if (Receivers->NumberOfBinsH != Antenna->GetNumberOfFrequencies())
  187. Receivers->SetNumberOfBinsH(Antenna->GetNumberOfFrequencies());
  188. break;
  189. case BOTH:
  190. if (Receivers->NumberOfBinsH != Antenna->GetNumberOfFrequencies())
  191. Receivers->SetNumberOfBinsH(Antenna->GetNumberOfFrequencies());
  192. if (Receivers->NumberOfBinsE != Antenna->GetNumberOfFrequencies())
  193. Receivers->SetNumberOfBinsE(Antenna->GetNumberOfFrequencies());
  194. break;
  195. }
  196. if (Antenna->GetName() == std::string("PolygonalWireAntenna") || Antenna->GetName() == std::string("TEMTransmitter") ) {
  197. icalc += 1;
  198. // Check to see if they are all on a plane? If so we can do this fast
  199. if (Antenna->IsHorizontallyPlanar() && HankelType == ANDERSON801) {
  200. for (int ifreq=0; ifreq<Antenna->GetNumberOfFrequencies();++ifreq) {
  201. Real wavef = 2.*PI* Antenna->GetFrequency(ifreq);
  202. #ifdef LEMMAUSEOMP
  203. #pragma omp parallel
  204. {
  205. #endif
  206. auto Hankel = FHTAnderson801::NewSP();
  207. #ifdef LEMMAUSEOMP
  208. #pragma omp for schedule(static, 1)
  209. #endif
  210. for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
  211. auto AntCopy = static_cast<PolygonalWireAntenna*>(Antenna.get())->ClonePA();
  212. SolveLaggedTxRxPair(irec, Hankel.get(), wavef, ifreq, AntCopy.get());
  213. }
  214. //Receivers->ClearFields();
  215. #ifdef LEMMAUSEOMP
  216. }
  217. #endif
  218. }
  219. } else
  220. if (Receivers->GetNumberOfPoints() > Antenna->GetNumberOfFrequencies()) {
  221. //std::cout << "freq parallel #1" << std::endl;
  222. //** Progress display bar for long calculations */
  223. #ifdef HAVE_BOOST_PROGRESS
  224. if (progressbar) {
  225. disp = new boost::progress_display( Receivers->GetNumberOfPoints()*Antenna->GetNumberOfFrequencies() );
  226. }
  227. #endif
  228. // parallelise across receivers
  229. #ifdef LEMMAUSEOMP
  230. #pragma omp parallel
  231. #endif
  232. { // OpenMP Parallel Block
  233. // Since these antennas change we need a local copy for each
  234. // thread.
  235. auto AntCopy = static_cast<PolygonalWireAntenna*>(Antenna.get())->ClonePA();
  236. std::shared_ptr<HankelTransform> Hankel;
  237. switch (HankelType) {
  238. case ANDERSON801:
  239. Hankel = FHTAnderson801::NewSP();
  240. break;
  241. case CHAVE:
  242. Hankel = GQChave::NewSP();
  243. break;
  244. case FHTKEY201:
  245. Hankel = FHTKey201::NewSP();
  246. break;
  247. case FHTKEY101:
  248. Hankel = FHTKey101::NewSP();
  249. break;
  250. case FHTKEY51:
  251. Hankel = FHTKey51::NewSP();
  252. break;
  253. case QWEKEY:
  254. Hankel = QWEKey::NewSP();
  255. break;
  256. default:
  257. std::cerr << "Hankel transform cannot be created\n";
  258. exit(EXIT_FAILURE);
  259. }
  260. //for (int irec=tid; irec<Receivers->GetNumberOfPoints(); irec+=nthreads) {
  261. #ifdef LEMMAUSEOMP
  262. #pragma omp for schedule(static, 1) //nowait
  263. #endif
  264. for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
  265. if (!Receivers->GetMask(irec)) {
  266. AntCopy->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
  267. for (int idip=0; idip<AntCopy->GetNumberOfDipoles(); ++idip) {
  268. auto tDipole = AntCopy->GetDipoleSource(idip);
  269. //#ifdef LEMMAUSEOMP
  270. //#pragma omp for schedule(static, 1)
  271. //#endif
  272. for (int ifreq=0; ifreq<tDipole->GetNumberOfFrequencies();
  273. ++ifreq) {
  274. // Propogation constant in free space
  275. Real wavef = tDipole->GetAngularFrequency(ifreq) *
  276. std::sqrt(MU0*EPSILON0);
  277. SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
  278. } // freq loop
  279. } // dipole loop
  280. } // mask
  281. //std::cout << "Normal Path\n";
  282. //std::cout << Receivers->GetHfield(0, irec) << std::endl;
  283. //if (irec == 1) exit(0);
  284. #ifdef HAVE_BOOST_PROGRESS
  285. if (progressbar) ++(*disp);
  286. #endif
  287. } // receiver loop
  288. } // OMP_PARALLEL BLOCK
  289. } else if (Antenna->GetNumberOfFrequencies() > 8) {
  290. // parallel across frequencies
  291. //std::cout << "freq parallel #2" << std::endl;
  292. for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
  293. if (!Receivers->GetMask(irec)) {
  294. static_cast<PolygonalWireAntenna*>(Antenna.get())->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
  295. #ifdef LEMMAUSEOMP
  296. #pragma omp parallel
  297. #endif
  298. { // OpenMP Parallel Block
  299. std::shared_ptr<HankelTransform> Hankel;
  300. switch (HankelType) {
  301. case ANDERSON801:
  302. Hankel = FHTAnderson801::NewSP();
  303. break;
  304. case CHAVE:
  305. Hankel = GQChave::NewSP();
  306. break;
  307. case FHTKEY201:
  308. Hankel = FHTKey201::NewSP();
  309. break;
  310. case FHTKEY101:
  311. Hankel = FHTKey101::NewSP();
  312. break;
  313. case FHTKEY51:
  314. Hankel = FHTKey51::NewSP();
  315. break;
  316. case QWEKEY:
  317. Hankel = QWEKey::NewSP();
  318. break;
  319. default:
  320. std::cerr << "Hankel transform cannot be created\n";
  321. exit(EXIT_FAILURE);
  322. }
  323. #ifdef LEMMAUSEOMP
  324. #pragma omp for schedule(static, 1)
  325. #endif
  326. for (int ifreq=0; ifreq<Antenna->GetNumberOfFrequencies(); ++ifreq) {
  327. for (int idip=0; idip<Antenna->GetNumberOfDipoles(); ++idip) {
  328. auto tDipole = Antenna->GetDipoleSource(idip);
  329. // Propogation constant in free space
  330. Real wavef = tDipole->GetAngularFrequency(ifreq) *
  331. std::sqrt(MU0*EPSILON0);
  332. SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
  333. } // dipole loop
  334. } // frequency loop
  335. } // OMP_PARALLEL BLOCK
  336. } // mask loop
  337. #ifdef HAVE_BOOST_PROGRESS
  338. //if (Receivers->GetNumberOfPoints() > 100) {
  339. // ++ disp;
  340. //}
  341. #endif
  342. } // receiver loop
  343. //std::cout << "End freq parallel " << std::endl;
  344. } // Frequency Parallel
  345. else {
  346. //std::cout << "parallel across #3 " << std::endl;
  347. for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
  348. if (!Receivers->GetMask(irec)) {
  349. static_cast<PolygonalWireAntenna*>(Antenna.get())->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
  350. // std::cout << "Not Masked " << std::endl;
  351. // std::cout << "n Freqs " << Antenna->GetNumberOfFrequencies() << std::endl;
  352. // std::cout << "n Dipoles " << Antenna->GetNumberOfDipoles() << std::endl;
  353. // if ( !Antenna->GetNumberOfDipoles() ) {
  354. // std::cout << "NO DIPOLES!!!!!!!!!!!!!!!!!!!!!!!!!!\n";
  355. // // std::cout << "rec location " << Receivers->GetLocation(irec) << std::endl;
  356. // // }
  357. #ifdef LEMMAUSEOMP
  358. #pragma omp parallel
  359. #endif
  360. { // OpenMP Parallel Block
  361. std::shared_ptr<HankelTransform> Hankel;
  362. switch (HankelType) {
  363. case ANDERSON801:
  364. Hankel = FHTAnderson801::NewSP();
  365. break;
  366. case CHAVE:
  367. Hankel = GQChave::NewSP();
  368. break;
  369. case FHTKEY201:
  370. Hankel = FHTKey201::NewSP();
  371. break;
  372. case FHTKEY101:
  373. Hankel = FHTKey101::NewSP();
  374. break;
  375. case FHTKEY51:
  376. Hankel = FHTKey51::NewSP();
  377. break;
  378. case QWEKEY:
  379. Hankel = QWEKey::NewSP();
  380. break;
  381. default:
  382. std::cerr << "Hankel transform cannot be created\n";
  383. exit(EXIT_FAILURE);
  384. }
  385. for (int ifreq=0; ifreq<Antenna->GetNumberOfFrequencies(); ++ifreq) {
  386. #ifdef LEMMAUSEOMP
  387. #pragma omp for schedule(static, 1)
  388. #endif
  389. for (int idip=0; idip<Antenna->GetNumberOfDipoles(); ++idip) {
  390. //#pragma omp critical
  391. //{
  392. //cout << "idip=" << idip << "\tthread num=" << omp_get_thread_num() << '\n';
  393. //}
  394. auto tDipole = Antenna->GetDipoleSource(idip);
  395. // Propogation constant in free space
  396. Real wavef = tDipole->GetAngularFrequency(ifreq) *
  397. std::sqrt(MU0*EPSILON0);
  398. SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
  399. } // dipole loop
  400. } // frequency loop
  401. } // OMP_PARALLEL BLOCK
  402. } // mask loop
  403. #ifdef HAVE_BOOST_PROGRESS
  404. //if (Receivers->GetNumberOfPoints() > 100) {
  405. // ++ disp;
  406. //}
  407. #endif
  408. } // receiver loop
  409. } // Polygonal parallel logic
  410. } else {
  411. std::cerr << "Lemma with WireAntenna class is currently broken"
  412. << " fix or use PolygonalWireAntenna\n" << std::endl;
  413. exit(EXIT_FAILURE);
  414. // TODO, getting wrong answer, curiously worKernel->GetKs() with MakeCalc, maybe
  415. // a threading issue, use SolveSingleTxRxPair maype instead of call
  416. // to MakeCalc3? !!!
  417. for (int idip=0; idip<Antenna->GetNumberOfDipoles(); ++idip) {
  418. this->Dipole = Antenna->GetDipoleSource(idip);
  419. MakeCalc3();
  420. //++disp;
  421. }
  422. this->Dipole = nullptr;
  423. }
  424. #ifdef HAVE_BOOST_PROGRESS
  425. if (progressbar) {
  426. delete disp;
  427. }
  428. #endif
  429. }
  430. #ifdef KIHALEE_EM1D
  431. void EMEarth1D::MakeCalc() {
  432. int itype; // 1 = elec, 2 = mag
  433. switch (this->Dipole->GetDipoleSourceType()) {
  434. case (GROUNDEDELECTRICDIPOLE) :
  435. itype = 1;
  436. break;
  437. case (MAGNETICDIPOLE) :
  438. itype = 2;
  439. break;
  440. case (UNGROUNDEDELECTRICDIPOLE) :
  441. std::cerr << "Fortran routine cannot calculate ungrounded"
  442. "electric dipole\n";
  443. default:
  444. throw NonValidDipoleType();
  445. }
  446. int ipol ;
  447. Vector3r Pol = this->Dipole->GetPolarisation();
  448. if (std::abs(Pol[0]-1) < 1e-5) {
  449. ipol = 1;
  450. } else if (std::abs(Pol[1]-1) < 1e-5) {
  451. ipol = 2;
  452. } else if (std::abs(Pol[2]-1) < 1e-5) {
  453. ipol = 3;
  454. } else {
  455. std::cerr << "Fortran routine cannot calculate arbitrary "
  456. "dipole polarisation, set to x, y, or z\n";
  457. }
  458. int nlay = Earth->GetNumberOfNonAirLayers();
  459. if (nlay > MAXLAYERS) {
  460. std::cerr << "FORTRAN CODE CAN ONLY HANDLE " << MAXLAYERS
  461. << " LAYERS\n";
  462. throw EarthModelWithMoreThanMaxLayers();
  463. }
  464. int nfreq = 1; // number of freqs
  465. int nfield; // field output 1 = elec, 2 = mag, 3 = both
  466. switch (FieldsToCalculate) {
  467. case E:
  468. nfield = 1;
  469. break;
  470. case H:
  471. nfield = 2;
  472. break;
  473. case BOTH:
  474. nfield = 3;
  475. break;
  476. default:
  477. throw 7;
  478. }
  479. int nres = Receivers->GetNumberOfPoints();
  480. int jtype = 3; // form ouf output,
  481. // 1 = horizontal,
  482. // 2 = down hole,
  483. // 3 = freq sounding
  484. // 4 = down hole logging
  485. int jgamma = 0; // Units 0 = MKS (H->A/m and E->V/m)
  486. // 1 = h->Gammas E->V/m
  487. double acc = 0.; // Tolerance
  488. // TODO, fix FORTRAN calls so these arrays can be nlay long, not
  489. // MAXLAYERS.
  490. // Model Parameters
  491. double *dep = new double[MAXLAYERS];
  492. dep[0] = 0.; // We always say air starts at 0
  493. for (int ilay=1; ilay<Earth->GetNumberOfLayers(); ++ilay) {
  494. dep[ilay] = dep[ilay-1] + Earth->GetLayerThickness(ilay);
  495. //std::cout << "Depth " << dep[ilay] << std::endl;
  496. }
  497. std::complex<double> *sig = new std::complex<double> [MAXLAYERS];
  498. for (int ilay=1; ilay<=nlay; ++ilay) {
  499. sig[ilay-1] = (std::complex<double>)(Earth->GetLayerConductivity(ilay));
  500. }
  501. // TODO, pass these into Fortran call, and return Cole-Cole model
  502. // parameters. Right now this does nothing
  503. //std::complex<double> *sus = new std::complex<double>[MAXLAYERS];
  504. //std::complex<double> *epr = new std::complex<double>[MAXLAYERS];
  505. // Cole-Cole model stuff
  506. double *susl = new double[MAXLAYERS];
  507. for (int ilay=1; ilay<=nlay; ++ilay) {
  508. susl[ilay-1] = Earth->GetLayerLowFreqSusceptibility(ilay);
  509. }
  510. double *sush = new double[MAXLAYERS];
  511. for (int ilay=1; ilay<=nlay; ++ilay) {
  512. sush[ilay-1] = Earth->GetLayerHighFreqSusceptibility(ilay);
  513. }
  514. double *sustau = new double[MAXLAYERS];
  515. for (int ilay=1; ilay<=nlay; ++ilay) {
  516. sustau[ilay-1] = Earth->GetLayerTauSusceptibility(ilay);
  517. }
  518. double *susalp = new double[MAXLAYERS];
  519. for (int ilay=1; ilay<=nlay; ++ilay) {
  520. susalp[ilay-1] = Earth->GetLayerBreathSusceptibility(ilay);
  521. }
  522. double *eprl = new double[MAXLAYERS];
  523. for (int ilay=1; ilay<=nlay; ++ilay) {
  524. eprl[ilay-1] = Earth->GetLayerLowFreqPermitivity(ilay);
  525. }
  526. double *eprh = new double[MAXLAYERS];
  527. for (int ilay=1; ilay<=nlay; ++ilay) {
  528. eprh[ilay-1] = Earth->GetLayerHighFreqPermitivity(ilay);
  529. }
  530. double *eprtau = new double[MAXLAYERS];
  531. for (int ilay=1; ilay<=nlay; ++ilay) {
  532. eprtau[ilay-1] = Earth->GetLayerTauPermitivity(ilay);
  533. }
  534. double *epralp = new double[MAXLAYERS];
  535. for (int ilay=1; ilay<=nlay; ++ilay) {
  536. epralp[ilay-1] = Earth->GetLayerBreathPermitivity(ilay);
  537. }
  538. // Freq stuff
  539. double finit = Dipole->GetFrequency(0); //(1000); // Starting freq
  540. double flimit = Dipole->GetFrequency(0); //(1000); // max freq
  541. double dlimit = Dipole->GetFrequency(0); //(1000); // difusion limit
  542. double lfinc(1); // no. freq per decade
  543. // tx location jtype != 4
  544. double txx = Dipole->GetLocation(0); // (0.);
  545. double txy = Dipole->GetLocation(1); // (0.);
  546. double txz = Dipole->GetLocation(2); // (0.);
  547. // rx position
  548. // TODO, fix Fortran program to not waste this memory
  549. // maybe
  550. const int MAXREC = 15;
  551. double *rxx = new double [MAXREC];
  552. double *rxy = new double [MAXREC];
  553. double *rxz = new double [MAXREC];
  554. std::complex<double> *ex = new std::complex<double>[MAXREC];
  555. std::complex<double> *ey = new std::complex<double>[MAXREC];
  556. std::complex<double> *ez = new std::complex<double>[MAXREC];
  557. std::complex<double> *hx = new std::complex<double>[MAXREC];
  558. std::complex<double> *hy = new std::complex<double>[MAXREC];
  559. std::complex<double> *hz = new std::complex<double>[MAXREC];
  560. int nres2 = MAXREC;
  561. int ii=0;
  562. for (ii=0; ii<nres-MAXREC; ii+=MAXREC) {
  563. for (int ir=0; ir<MAXREC; ++ir) {
  564. //Vector3r pos = Receivers->GetLocation(ii+ir);
  565. rxx[ir] = Receivers->GetLocation(ii+ir)[0];
  566. rxy[ir] = Receivers->GetLocation(ii+ir)[1];
  567. rxz[ir] = Receivers->GetLocation(ii+ir)[2];
  568. }
  569. em1dcall_(itype, ipol, nlay, nfreq, nfield, nres2, jtype,
  570. jgamma, acc, dep, sig, susl, sush, sustau, susalp,
  571. eprl, eprh, eprtau, epralp, finit, flimit, dlimit,
  572. lfinc, txx, txy, txz, rxx, rxy, rxz, ex, ey, ez,
  573. hx, hy, hz);
  574. // Scale By Moment
  575. for (int ir=0; ir<MAXREC; ++ir) {
  576. ex[ir] *= Dipole->GetMoment();
  577. ey[ir] *= Dipole->GetMoment();
  578. ez[ir] *= Dipole->GetMoment();
  579. hx[ir] *= Dipole->GetMoment();
  580. hy[ir] *= Dipole->GetMoment();
  581. hz[ir] *= Dipole->GetMoment();
  582. // Append values instead of setting them
  583. this->Receivers->AppendEfield(0, ii+ir, (Complex)(ex[ir]),
  584. (Complex)(ey[ir]),
  585. (Complex)(ez[ir]) );
  586. this->Receivers->AppendHfield(0, ii+ir, (Complex)(hx[ir]),
  587. (Complex)(hy[ir]),
  588. (Complex)(hz[ir]) );
  589. }
  590. }
  591. //ii += MAXREC;
  592. nres2 = 0;
  593. // Perform last positions
  594. for (int ir=0; ir<nres-ii; ++ir) {
  595. rxx[ir] = Receivers->GetLocation(ii+ir)[0];
  596. rxy[ir] = Receivers->GetLocation(ii+ir)[1];
  597. rxz[ir] = Receivers->GetLocation(ii+ir)[2];
  598. ++nres2;
  599. }
  600. em1dcall_(itype, ipol, nlay, nfreq, nfield, nres2, jtype,
  601. jgamma, acc, dep, sig, susl, sush, sustau, susalp,
  602. eprl, eprh, eprtau, epralp, finit, flimit, dlimit,
  603. lfinc, txx, txy, txz, rxx, rxy, rxz, ex, ey, ez,
  604. hx, hy, hz);
  605. // Scale By Moment
  606. for (int ir=0; ir<nres-ii; ++ir) {
  607. ex[ir] *= Dipole->GetMoment();
  608. ey[ir] *= Dipole->GetMoment();
  609. ez[ir] *= Dipole->GetMoment();
  610. hx[ir] *= Dipole->GetMoment();
  611. hy[ir] *= Dipole->GetMoment();
  612. hz[ir] *= Dipole->GetMoment();
  613. // Append values instead of setting them
  614. this->Receivers->AppendEfield(0, ii+ir, (Complex)(ex[ir]),
  615. (Complex)(ey[ir]),
  616. (Complex)(ez[ir]) );
  617. this->Receivers->AppendHfield(0, ii+ir, (Complex)(hx[ir]),
  618. (Complex)(hy[ir]),
  619. (Complex)(hz[ir]) );
  620. }
  621. delete [] sig;
  622. delete [] dep;
  623. //delete [] sus;
  624. //delete [] epr;
  625. delete [] susl;
  626. delete [] sush;
  627. delete [] susalp;
  628. delete [] sustau;
  629. delete [] eprl;
  630. delete [] eprh;
  631. delete [] epralp;
  632. delete [] eprtau;
  633. delete [] rxx;
  634. delete [] rxy;
  635. delete [] rxz;
  636. delete [] ex;
  637. delete [] ey;
  638. delete [] ez;
  639. delete [] hx;
  640. delete [] hy;
  641. delete [] hz;
  642. }
  643. #endif
  644. void EMEarth1D::SolveSingleTxRxPair (const int &irec, HankelTransform *Hankel, const Real &wavef, const int &ifreq,
  645. DipoleSource *tDipole) {
  646. ++icalcinner;
  647. Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
  648. tDipole->SetKernels(ifreq, FieldsToCalculate, Receivers, irec, Earth);
  649. Hankel->ComputeRelated( rho, tDipole->GetKernelManager() );
  650. tDipole->UpdateFields( ifreq, Hankel, wavef );
  651. }
  652. // void EMEarth1D::SolveSingleTxRxPair (const int &irec, std::shared_ptr<HankelTransform> Hankel, const Real &wavef, const int &ifreq,
  653. // std::shared_ptr<DipoleSource> tDipole) {
  654. // ++icalcinner;
  655. // Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
  656. // tDipole->SetKernels(ifreq, FieldsToCalculate, Receivers, irec, Earth);
  657. // //Hankel->ComputeRelated( rho, tDipole->GetKernelManager() );
  658. // //tDipole->UpdateFields( ifreq, Hankel, wavef );
  659. // }
  660. void EMEarth1D::SolveLaggedTxRxPair(const int &irec, FHTAnderson801* Hankel,
  661. const Real &wavef, const int &ifreq, PolygonalWireAntenna* antenna) {
  662. antenna->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
  663. // Determine the min and max arguments
  664. Real rhomin = 1e9;
  665. Real rhomax = 1e-9;
  666. for (int idip=0; idip<antenna->GetNumberOfDipoles(); ++idip) {
  667. auto tDipole = antenna->GetDipoleSource(idip);
  668. Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
  669. rhomin = std::min(rhomin, rho);
  670. rhomax = std::max(rhomax, rho);
  671. }
  672. //std::cout << "rhomin\t" << rhomin << "\trhomax" << rhomax << std::endl;
  673. // Determine number of lagged convolutions to do
  674. // TODO, can Hankel2 adjust the lagg spacing safely?
  675. int nlag = 1; // We need an extra for some reason for stability
  676. Real lrho ( 1.01* rhomax );
  677. while ( lrho > rhomin ) {
  678. nlag += 1;
  679. lrho *= Hankel->GetABSER();
  680. }
  681. //int nlag = rhomin
  682. auto tDipole = antenna->GetDipoleSource(0);
  683. tDipole->SetKernels(ifreq, FieldsToCalculate, Receivers, irec, Earth);
  684. // Instead we should pass the antenna into this so that Hankel hass all the rho arguments...
  685. Hankel->ComputeLaggedRelated( 1.01* rhomax, nlag, tDipole->GetKernelManager() );
  686. //std::cout << Hankel->GetAnswer() << std::endl;
  687. //std::cout << Hankel->GetArg() << std::endl;
  688. // Sort the dipoles by rho
  689. for (int idip=0; idip<antenna->GetNumberOfDipoles(); ++idip) {
  690. //for (int idip=0; idip<1; ++idip) {
  691. auto tDipole = antenna->GetDipoleSource(idip);
  692. tDipole->SetKernels(ifreq, FieldsToCalculate, Receivers, irec, Earth);
  693. // Pass Hankel2 a message here so it knows which one to return in Zgauss!
  694. Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
  695. //std::cout << " in Lagged " << rho << "\t" << rhomin << "\t" << rhomax << std::endl;
  696. Hankel->SetLaggedArg( rho );
  697. //std::cout << "out Lagged" << std::endl;
  698. tDipole->UpdateFields( ifreq, Hankel, wavef );
  699. }
  700. //std::cout << "Spline\n";
  701. //std::cout << Receivers->GetHfield(0, irec) << std::endl;
  702. }
  703. //////////////////////////////////////////////////////////
  704. // Thread safe OO Reimplimentation of KiHand's
  705. // EM1DNEW.for programme
  706. void EMEarth1D::MakeCalc3() {
  707. if ( Dipole == nullptr ) throw NullDipoleSource();
  708. if (Earth == nullptr) throw NullEarth();
  709. if (Receivers == nullptr) throw NullReceivers();
  710. #ifdef LEMMAUSEOMP
  711. #pragma omp parallel
  712. #endif
  713. { // OpenMP Parallel Block
  714. #ifdef LEMMAUSEOMP
  715. int tid = omp_get_thread_num();
  716. int nthreads = omp_get_num_threads();
  717. #else
  718. int tid=0;
  719. int nthreads=1;
  720. #endif
  721. auto tDipole = Dipole->Clone();
  722. std::shared_ptr<HankelTransform> Hankel;
  723. switch (HankelType) {
  724. case ANDERSON801:
  725. Hankel = FHTAnderson801::NewSP();
  726. break;
  727. case CHAVE:
  728. Hankel = GQChave::NewSP();
  729. break;
  730. case FHTKEY201:
  731. Hankel = FHTKey201::NewSP();
  732. break;
  733. case FHTKEY101:
  734. Hankel = FHTKey101::NewSP();
  735. break;
  736. case FHTKEY51:
  737. Hankel = FHTKey51::NewSP();
  738. break;
  739. case QWEKEY:
  740. Hankel = QWEKey::NewSP();
  741. break;
  742. default:
  743. std::cerr << "Hankel transform cannot be created\n";
  744. exit(EXIT_FAILURE);
  745. }
  746. if ( tDipole->GetNumberOfFrequencies() < Receivers->GetNumberOfPoints() ) {
  747. for (int ifreq=0; ifreq<tDipole->GetNumberOfFrequencies(); ++ifreq) {
  748. // Propogation constant in free space being input to Hankel
  749. Real wavef = tDipole->GetAngularFrequency(ifreq) * std::sqrt(MU0*EPSILON0);
  750. for (int irec=tid; irec<Receivers->GetNumberOfPoints(); irec+=nthreads) {
  751. SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
  752. }
  753. }
  754. } else {
  755. for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
  756. for (int ifreq=tid; ifreq<tDipole->GetNumberOfFrequencies(); ifreq+=nthreads) {
  757. // Propogation constant in free space being input to Hankel
  758. Real wavef = tDipole->GetAngularFrequency(ifreq) * std::sqrt(MU0*EPSILON0);
  759. SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
  760. }
  761. }
  762. }
  763. } // OpenMP Parallel Block
  764. }
  765. NullReceivers::NullReceivers() :
  766. runtime_error("nullptr RECEIVERS") {}
  767. NullAntenna::NullAntenna() :
  768. runtime_error("nullptr ANTENNA") {}
  769. NullInstrument::NullInstrument(LemmaObject* ptr) :
  770. runtime_error("nullptr INSTRUMENT") {
  771. std::cout << "Thrown by instance of "
  772. << ptr->GetName() << std::endl;
  773. }
  774. DipoleSourceSpecifiedForWireAntennaCalc::
  775. DipoleSourceSpecifiedForWireAntennaCalc() :
  776. runtime_error("DIPOLE SOURCE SPECIFIED FOR WIRE ANTENNA CALC"){}
  777. } // end of Lemma Namespace