Lemma is an Electromagnetics API
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

EMEarth1D.cpp 37KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900
  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 HAVEBOOSTPROGRESS
  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 HAVEBOOSTPROGRESS
  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. /* TODO FIX THIS ISSUES */
  200. if (Antenna->IsHorizontallyPlanar() && HankelType == ANDERSON801) {
  201. //std::cout << "Lag baby lag" << std::endl;
  202. for (int ifreq=0; ifreq<Antenna->GetNumberOfFrequencies();++ifreq) {
  203. //std::cout << "Num Recs" << Receivers->GetNumberOfPoints() << std::endl;
  204. Real wavef = 2.*PI* Antenna->GetFrequency(ifreq);
  205. #ifdef LEMMAUSEOMP
  206. #pragma omp parallel
  207. {
  208. #endif
  209. auto Hankel = FHTAnderson801::NewSP();
  210. #ifdef LEMMAUSEOMP
  211. #pragma omp for schedule(static, 1)
  212. #endif
  213. for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
  214. //for (int irec=0; irec<2; ++irec) { // TODO FIXME BELO
  215. auto AntCopy = static_cast<PolygonalWireAntenna*>(Antenna.get())->ClonePA();
  216. SolveLaggedTxRxPair(irec, Hankel.get(), wavef, ifreq, AntCopy.get());
  217. //exit(0);
  218. }
  219. //Receivers->ClearFields(); // FIXME DEBUG TODO
  220. #ifdef LEMMAUSEOMP
  221. }
  222. #endif
  223. }
  224. } else
  225. if (Receivers->GetNumberOfPoints() > Antenna->GetNumberOfFrequencies()) {
  226. //std::cout << "freq parallel #1" << std::endl;
  227. //** Progress display bar for long calculations */
  228. #ifdef HAVEBOOSTPROGRESS
  229. if (progressbar) {
  230. disp = new boost::progress_display( Receivers->GetNumberOfPoints()*Antenna->GetNumberOfFrequencies() );
  231. }
  232. #endif
  233. // parallelise across receivers
  234. #ifdef LEMMAUSEOMP
  235. #pragma omp parallel
  236. #endif
  237. { // OpenMP Parallel Block
  238. // Since these antennas change we need a local copy for each
  239. // thread.
  240. auto AntCopy = static_cast<PolygonalWireAntenna*>(Antenna.get())->ClonePA();
  241. std::shared_ptr<HankelTransform> Hankel;
  242. switch (HankelType) {
  243. case ANDERSON801:
  244. Hankel = FHTAnderson801::NewSP();
  245. break;
  246. case CHAVE:
  247. Hankel = GQChave::NewSP();
  248. break;
  249. case FHTKEY201:
  250. Hankel = FHTKey201::NewSP();
  251. break;
  252. case FHTKEY101:
  253. Hankel = FHTKey101::NewSP();
  254. break;
  255. case FHTKEY51:
  256. Hankel = FHTKey51::NewSP();
  257. break;
  258. case QWEKEY:
  259. Hankel = QWEKey::NewSP();
  260. break;
  261. default:
  262. std::cerr << "Hankel transform cannot be created\n";
  263. exit(EXIT_FAILURE);
  264. }
  265. //for (int irec=tid; irec<Receivers->GetNumberOfPoints(); irec+=nthreads) {
  266. #ifdef LEMMAUSEOMP
  267. #pragma omp for schedule(static, 1) //nowait
  268. #endif
  269. for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
  270. if (!Receivers->GetMask(irec)) {
  271. AntCopy->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
  272. for (int idip=0; idip<AntCopy->GetNumberOfDipoles(); ++idip) {
  273. auto tDipole = AntCopy->GetDipoleSource(idip);
  274. //#ifdef LEMMAUSEOMP
  275. //#pragma omp for schedule(static, 1)
  276. //#endif
  277. for (int ifreq=0; ifreq<tDipole->GetNumberOfFrequencies();
  278. ++ifreq) {
  279. // Propogation constant in free space
  280. Real wavef = tDipole->GetAngularFrequency(ifreq) *
  281. std::sqrt(MU0*EPSILON0);
  282. SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
  283. } // freq loop
  284. } // dipole loop
  285. } // mask
  286. //std::cout << "Normal Path\n";
  287. //std::cout << Receivers->GetHfield(0, irec) << std::endl;
  288. //if (irec == 1) exit(0);
  289. #ifdef HAVEBOOSTPROGRESS
  290. if (progressbar) ++(*disp);
  291. #endif
  292. } // receiver loop
  293. } // OMP_PARALLEL BLOCK
  294. } else if (Antenna->GetNumberOfFrequencies() > 8) {
  295. // parallel across frequencies
  296. //std::cout << "freq parallel #2" << std::endl;
  297. for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
  298. if (!Receivers->GetMask(irec)) {
  299. static_cast<PolygonalWireAntenna*>(Antenna.get())->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
  300. #ifdef LEMMAUSEOMP
  301. #pragma omp parallel
  302. #endif
  303. { // OpenMP Parallel Block
  304. std::shared_ptr<HankelTransform> Hankel;
  305. switch (HankelType) {
  306. case ANDERSON801:
  307. Hankel = FHTAnderson801::NewSP();
  308. break;
  309. case CHAVE:
  310. Hankel = GQChave::NewSP();
  311. break;
  312. case FHTKEY201:
  313. Hankel = FHTKey201::NewSP();
  314. break;
  315. case FHTKEY101:
  316. Hankel = FHTKey101::NewSP();
  317. break;
  318. case FHTKEY51:
  319. Hankel = FHTKey51::NewSP();
  320. break;
  321. case QWEKEY:
  322. Hankel = QWEKey::NewSP();
  323. break;
  324. default:
  325. std::cerr << "Hankel transform cannot be created\n";
  326. exit(EXIT_FAILURE);
  327. }
  328. #ifdef LEMMAUSEOMP
  329. #pragma omp for schedule(static, 1)
  330. #endif
  331. for (int ifreq=0; ifreq<Antenna->GetNumberOfFrequencies(); ++ifreq) {
  332. for (int idip=0; idip<Antenna->GetNumberOfDipoles(); ++idip) {
  333. auto tDipole = Antenna->GetDipoleSource(idip);
  334. // Propogation constant in free space
  335. Real wavef = tDipole->GetAngularFrequency(ifreq) *
  336. std::sqrt(MU0*EPSILON0);
  337. SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
  338. } // dipole loop
  339. } // frequency loop
  340. } // OMP_PARALLEL BLOCK
  341. } // mask loop
  342. #ifdef HAVEBOOSTPROGRESS
  343. //if (Receivers->GetNumberOfPoints() > 100) {
  344. // ++ disp;
  345. //}
  346. #endif
  347. } // receiver loop
  348. //std::cout << "End freq parallel " << std::endl;
  349. } // Frequency Parallel
  350. else {
  351. //std::cout << "parallel across #3 " << std::endl;
  352. for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
  353. if (!Receivers->GetMask(irec)) {
  354. static_cast<PolygonalWireAntenna*>(Antenna.get())->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
  355. // std::cout << "Not Masked " << std::endl;
  356. // std::cout << "n Freqs " << Antenna->GetNumberOfFrequencies() << std::endl;
  357. // std::cout << "n Dipoles " << Antenna->GetNumberOfDipoles() << std::endl;
  358. // if ( !Antenna->GetNumberOfDipoles() ) {
  359. // std::cout << "NO DIPOLES!!!!!!!!!!!!!!!!!!!!!!!!!!\n";
  360. // // std::cout << "rec location " << Receivers->GetLocation(irec) << std::endl;
  361. // // }
  362. #ifdef LEMMAUSEOMP
  363. #pragma omp parallel
  364. #endif
  365. { // OpenMP Parallel Block
  366. std::shared_ptr<HankelTransform> Hankel;
  367. switch (HankelType) {
  368. case ANDERSON801:
  369. Hankel = FHTAnderson801::NewSP();
  370. break;
  371. case CHAVE:
  372. Hankel = GQChave::NewSP();
  373. break;
  374. case FHTKEY201:
  375. Hankel = FHTKey201::NewSP();
  376. break;
  377. case FHTKEY101:
  378. Hankel = FHTKey101::NewSP();
  379. break;
  380. case FHTKEY51:
  381. Hankel = FHTKey51::NewSP();
  382. break;
  383. case QWEKEY:
  384. Hankel = QWEKey::NewSP();
  385. break;
  386. default:
  387. std::cerr << "Hankel transform cannot be created\n";
  388. exit(EXIT_FAILURE);
  389. }
  390. for (int ifreq=0; ifreq<Antenna->GetNumberOfFrequencies(); ++ifreq) {
  391. #ifdef LEMMAUSEOMP
  392. #pragma omp for schedule(static, 1)
  393. #endif
  394. for (int idip=0; idip<Antenna->GetNumberOfDipoles(); ++idip) {
  395. //#pragma omp critical
  396. //{
  397. //cout << "idip=" << idip << "\tthread num=" << omp_get_thread_num() << '\n';
  398. //}
  399. auto tDipole = Antenna->GetDipoleSource(idip);
  400. // Propogation constant in free space
  401. Real wavef = tDipole->GetAngularFrequency(ifreq) *
  402. std::sqrt(MU0*EPSILON0);
  403. SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
  404. } // dipole loop
  405. } // frequency loop
  406. } // OMP_PARALLEL BLOCK
  407. } // mask loop
  408. #ifdef HAVEBOOSTPROGRESS
  409. //if (Receivers->GetNumberOfPoints() > 100) {
  410. // ++ disp;
  411. //}
  412. #endif
  413. } // receiver loop
  414. } // Polygonal parallel logic
  415. } else {
  416. std::cerr << "Lemma with WireAntenna class is currently broken"
  417. << " fix or use PolygonalWireAntenna\n" << std::endl;
  418. exit(EXIT_FAILURE);
  419. // TODO, getting wrong answer, curiously worKernel->GetKs() with MakeCalc, maybe
  420. // a threading issue, use SolveSingleTxRxPair maype instead of call
  421. // to MakeCalc3? !!!
  422. for (int idip=0; idip<Antenna->GetNumberOfDipoles(); ++idip) {
  423. this->Dipole = Antenna->GetDipoleSource(idip);
  424. MakeCalc3();
  425. //++disp;
  426. }
  427. this->Dipole = nullptr;
  428. }
  429. #ifdef HAVEBOOSTPROGRESS
  430. if (progressbar) {
  431. delete disp;
  432. }
  433. #endif
  434. }
  435. #ifdef KIHALEE_EM1D
  436. void EMEarth1D::MakeCalc() {
  437. int itype; // 1 = elec, 2 = mag
  438. switch (this->Dipole->GetDipoleSourceType()) {
  439. case (GROUNDEDELECTRICDIPOLE) :
  440. itype = 1;
  441. break;
  442. case (MAGNETICDIPOLE) :
  443. itype = 2;
  444. break;
  445. case (UNGROUNDEDELECTRICDIPOLE) :
  446. std::cerr << "Fortran routine cannot calculate ungrounded"
  447. "electric dipole\n";
  448. default:
  449. throw NonValidDipoleType();
  450. }
  451. int ipol ;
  452. Vector3r Pol = this->Dipole->GetPolarisation();
  453. if (std::abs(Pol[0]-1) < 1e-5) {
  454. ipol = 1;
  455. } else if (std::abs(Pol[1]-1) < 1e-5) {
  456. ipol = 2;
  457. } else if (std::abs(Pol[2]-1) < 1e-5) {
  458. ipol = 3;
  459. } else {
  460. std::cerr << "Fortran routine cannot calculate arbitrary "
  461. "dipole polarisation, set to x, y, or z\n";
  462. }
  463. int nlay = Earth->GetNumberOfNonAirLayers();
  464. if (nlay > MAXLAYERS) {
  465. std::cerr << "FORTRAN CODE CAN ONLY HANDLE " << MAXLAYERS
  466. << " LAYERS\n";
  467. throw EarthModelWithMoreThanMaxLayers();
  468. }
  469. int nfreq = 1; // number of freqs
  470. int nfield; // field output 1 = elec, 2 = mag, 3 = both
  471. switch (FieldsToCalculate) {
  472. case E:
  473. nfield = 1;
  474. break;
  475. case H:
  476. nfield = 2;
  477. break;
  478. case BOTH:
  479. nfield = 3;
  480. break;
  481. default:
  482. throw 7;
  483. }
  484. int nres = Receivers->GetNumberOfPoints();
  485. int jtype = 3; // form ouf output,
  486. // 1 = horizontal,
  487. // 2 = down hole,
  488. // 3 = freq sounding
  489. // 4 = down hole logging
  490. int jgamma = 0; // Units 0 = MKS (H->A/m and E->V/m)
  491. // 1 = h->Gammas E->V/m
  492. double acc = 0.; // Tolerance
  493. // TODO, fix FORTRAN calls so these arrays can be nlay long, not
  494. // MAXLAYERS.
  495. // Model Parameters
  496. double *dep = new double[MAXLAYERS];
  497. dep[0] = 0.; // We always say air starts at 0
  498. for (int ilay=1; ilay<Earth->GetNumberOfLayers(); ++ilay) {
  499. dep[ilay] = dep[ilay-1] + Earth->GetLayerThickness(ilay);
  500. //std::cout << "Depth " << dep[ilay] << std::endl;
  501. }
  502. std::complex<double> *sig = new std::complex<double> [MAXLAYERS];
  503. for (int ilay=1; ilay<=nlay; ++ilay) {
  504. sig[ilay-1] = (std::complex<double>)(Earth->GetLayerConductivity(ilay));
  505. }
  506. // TODO, pass these into Fortran call, and return Cole-Cole model
  507. // parameters. Right now this does nothing
  508. //std::complex<double> *sus = new std::complex<double>[MAXLAYERS];
  509. //std::complex<double> *epr = new std::complex<double>[MAXLAYERS];
  510. // Cole-Cole model stuff
  511. double *susl = new double[MAXLAYERS];
  512. for (int ilay=1; ilay<=nlay; ++ilay) {
  513. susl[ilay-1] = Earth->GetLayerLowFreqSusceptibility(ilay);
  514. }
  515. double *sush = new double[MAXLAYERS];
  516. for (int ilay=1; ilay<=nlay; ++ilay) {
  517. sush[ilay-1] = Earth->GetLayerHighFreqSusceptibility(ilay);
  518. }
  519. double *sustau = new double[MAXLAYERS];
  520. for (int ilay=1; ilay<=nlay; ++ilay) {
  521. sustau[ilay-1] = Earth->GetLayerTauSusceptibility(ilay);
  522. }
  523. double *susalp = new double[MAXLAYERS];
  524. for (int ilay=1; ilay<=nlay; ++ilay) {
  525. susalp[ilay-1] = Earth->GetLayerBreathSusceptibility(ilay);
  526. }
  527. double *eprl = new double[MAXLAYERS];
  528. for (int ilay=1; ilay<=nlay; ++ilay) {
  529. eprl[ilay-1] = Earth->GetLayerLowFreqPermitivity(ilay);
  530. }
  531. double *eprh = new double[MAXLAYERS];
  532. for (int ilay=1; ilay<=nlay; ++ilay) {
  533. eprh[ilay-1] = Earth->GetLayerHighFreqPermitivity(ilay);
  534. }
  535. double *eprtau = new double[MAXLAYERS];
  536. for (int ilay=1; ilay<=nlay; ++ilay) {
  537. eprtau[ilay-1] = Earth->GetLayerTauPermitivity(ilay);
  538. }
  539. double *epralp = new double[MAXLAYERS];
  540. for (int ilay=1; ilay<=nlay; ++ilay) {
  541. epralp[ilay-1] = Earth->GetLayerBreathPermitivity(ilay);
  542. }
  543. // Freq stuff
  544. double finit = Dipole->GetFrequency(0); //(1000); // Starting freq
  545. double flimit = Dipole->GetFrequency(0); //(1000); // max freq
  546. double dlimit = Dipole->GetFrequency(0); //(1000); // difusion limit
  547. double lfinc(1); // no. freq per decade
  548. // tx location jtype != 4
  549. double txx = Dipole->GetLocation(0); // (0.);
  550. double txy = Dipole->GetLocation(1); // (0.);
  551. double txz = Dipole->GetLocation(2); // (0.);
  552. // rx position
  553. // TODO, fix Fortran program to not waste this memory
  554. // maybe
  555. const int MAXREC = 15;
  556. double *rxx = new double [MAXREC];
  557. double *rxy = new double [MAXREC];
  558. double *rxz = new double [MAXREC];
  559. std::complex<double> *ex = new std::complex<double>[MAXREC];
  560. std::complex<double> *ey = new std::complex<double>[MAXREC];
  561. std::complex<double> *ez = new std::complex<double>[MAXREC];
  562. std::complex<double> *hx = new std::complex<double>[MAXREC];
  563. std::complex<double> *hy = new std::complex<double>[MAXREC];
  564. std::complex<double> *hz = new std::complex<double>[MAXREC];
  565. int nres2 = MAXREC;
  566. int ii=0;
  567. for (ii=0; ii<nres-MAXREC; ii+=MAXREC) {
  568. for (int ir=0; ir<MAXREC; ++ir) {
  569. //Vector3r pos = Receivers->GetLocation(ii+ir);
  570. rxx[ir] = Receivers->GetLocation(ii+ir)[0];
  571. rxy[ir] = Receivers->GetLocation(ii+ir)[1];
  572. rxz[ir] = Receivers->GetLocation(ii+ir)[2];
  573. }
  574. em1dcall_(itype, ipol, nlay, nfreq, nfield, nres2, jtype,
  575. jgamma, acc, dep, sig, susl, sush, sustau, susalp,
  576. eprl, eprh, eprtau, epralp, finit, flimit, dlimit,
  577. lfinc, txx, txy, txz, rxx, rxy, rxz, ex, ey, ez,
  578. hx, hy, hz);
  579. // Scale By Moment
  580. for (int ir=0; ir<MAXREC; ++ir) {
  581. ex[ir] *= Dipole->GetMoment();
  582. ey[ir] *= Dipole->GetMoment();
  583. ez[ir] *= Dipole->GetMoment();
  584. hx[ir] *= Dipole->GetMoment();
  585. hy[ir] *= Dipole->GetMoment();
  586. hz[ir] *= Dipole->GetMoment();
  587. // Append values instead of setting them
  588. this->Receivers->AppendEfield(0, ii+ir, (Complex)(ex[ir]),
  589. (Complex)(ey[ir]),
  590. (Complex)(ez[ir]) );
  591. this->Receivers->AppendHfield(0, ii+ir, (Complex)(hx[ir]),
  592. (Complex)(hy[ir]),
  593. (Complex)(hz[ir]) );
  594. }
  595. }
  596. //ii += MAXREC;
  597. nres2 = 0;
  598. // Perform last positions
  599. for (int ir=0; ir<nres-ii; ++ir) {
  600. rxx[ir] = Receivers->GetLocation(ii+ir)[0];
  601. rxy[ir] = Receivers->GetLocation(ii+ir)[1];
  602. rxz[ir] = Receivers->GetLocation(ii+ir)[2];
  603. ++nres2;
  604. }
  605. em1dcall_(itype, ipol, nlay, nfreq, nfield, nres2, jtype,
  606. jgamma, acc, dep, sig, susl, sush, sustau, susalp,
  607. eprl, eprh, eprtau, epralp, finit, flimit, dlimit,
  608. lfinc, txx, txy, txz, rxx, rxy, rxz, ex, ey, ez,
  609. hx, hy, hz);
  610. // Scale By Moment
  611. for (int ir=0; ir<nres-ii; ++ir) {
  612. ex[ir] *= Dipole->GetMoment();
  613. ey[ir] *= Dipole->GetMoment();
  614. ez[ir] *= Dipole->GetMoment();
  615. hx[ir] *= Dipole->GetMoment();
  616. hy[ir] *= Dipole->GetMoment();
  617. hz[ir] *= Dipole->GetMoment();
  618. // Append values instead of setting them
  619. this->Receivers->AppendEfield(0, ii+ir, (Complex)(ex[ir]),
  620. (Complex)(ey[ir]),
  621. (Complex)(ez[ir]) );
  622. this->Receivers->AppendHfield(0, ii+ir, (Complex)(hx[ir]),
  623. (Complex)(hy[ir]),
  624. (Complex)(hz[ir]) );
  625. }
  626. delete [] sig;
  627. delete [] dep;
  628. //delete [] sus;
  629. //delete [] epr;
  630. delete [] susl;
  631. delete [] sush;
  632. delete [] susalp;
  633. delete [] sustau;
  634. delete [] eprl;
  635. delete [] eprh;
  636. delete [] epralp;
  637. delete [] eprtau;
  638. delete [] rxx;
  639. delete [] rxy;
  640. delete [] rxz;
  641. delete [] ex;
  642. delete [] ey;
  643. delete [] ez;
  644. delete [] hx;
  645. delete [] hy;
  646. delete [] hz;
  647. }
  648. #endif
  649. void EMEarth1D::SolveSingleTxRxPair (const int &irec, HankelTransform *Hankel, const Real &wavef, const int &ifreq,
  650. DipoleSource *tDipole) {
  651. ++icalcinner;
  652. Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
  653. tDipole->SetKernels(ifreq, FieldsToCalculate, Receivers, irec, Earth);
  654. Hankel->ComputeRelated( rho, tDipole->GetKernelManager() );
  655. tDipole->UpdateFields( ifreq, Hankel, wavef );
  656. }
  657. // void EMEarth1D::SolveSingleTxRxPair (const int &irec, std::shared_ptr<HankelTransform> Hankel, const Real &wavef, const int &ifreq,
  658. // std::shared_ptr<DipoleSource> tDipole) {
  659. // ++icalcinner;
  660. // Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
  661. // tDipole->SetKernels(ifreq, FieldsToCalculate, Receivers, irec, Earth);
  662. // //Hankel->ComputeRelated( rho, tDipole->GetKernelManager() );
  663. // //tDipole->UpdateFields( ifreq, Hankel, wavef );
  664. // }
  665. void EMEarth1D::SolveLaggedTxRxPair(const int &irec, FHTAnderson801* Hankel,
  666. const Real &wavef, const int &ifreq, PolygonalWireAntenna* antenna) {
  667. antenna->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
  668. // Determine the min and max arguments
  669. Real rhomin = 1e9;
  670. Real rhomax = 1e-9;
  671. for (int idip=0; idip<antenna->GetNumberOfDipoles(); ++idip) {
  672. auto tDipole = antenna->GetDipoleSource(idip);
  673. Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
  674. rhomin = std::min(rhomin, rho);
  675. rhomax = std::max(rhomax, rho);
  676. }
  677. //std::cout << "rhomin\t" << rhomin << "\trhomax" << rhomax << std::endl;
  678. // Determine number of lagged convolutions to do
  679. // TODO, can Hankel2 adjust the lagg spacing safely?
  680. int nlag = 1; // We need an extra for some reason for stability
  681. Real lrho ( 1.01* rhomax );
  682. while ( lrho > rhomin ) {
  683. nlag += 1;
  684. lrho *= Hankel->GetABSER();
  685. }
  686. //int nlag = rhomin
  687. auto tDipole = antenna->GetDipoleSource(0);
  688. tDipole->SetKernels(ifreq, FieldsToCalculate, Receivers, irec, Earth);
  689. // Instead we should pass the antenna into this so that Hankel hass all the rho arguments...
  690. Hankel->ComputeLaggedRelated( 1.01* rhomax, nlag, tDipole->GetKernelManager() );
  691. //std::cout << Hankel->GetAnswer() << std::endl;
  692. //std::cout << Hankel->GetArg() << std::endl;
  693. // Sort the dipoles by rho
  694. for (int idip=0; idip<antenna->GetNumberOfDipoles(); ++idip) {
  695. //for (int idip=0; idip<1; ++idip) {
  696. auto tDipole = antenna->GetDipoleSource(idip);
  697. tDipole->SetKernels(ifreq, FieldsToCalculate, Receivers, irec, Earth);
  698. // Pass Hankel2 a message here so it knows which one to return in Zgauss!
  699. Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
  700. //std::cout << " in Lagged " << rho << "\t" << rhomin << "\t" << rhomax << std::endl;
  701. Hankel->SetLaggedArg( rho );
  702. //std::cout << "out Lagged" << std::endl;
  703. tDipole->UpdateFields( ifreq, Hankel, wavef );
  704. }
  705. //std::cout << "Spline\n";
  706. //std::cout << Receivers->GetHfield(0, irec) << std::endl;
  707. }
  708. //////////////////////////////////////////////////////////
  709. // Thread safe OO Reimplimentation of KiHand's
  710. // EM1DNEW.for programme
  711. void EMEarth1D::MakeCalc3() {
  712. if ( Dipole == nullptr ) throw NullDipoleSource();
  713. if (Earth == nullptr) throw NullEarth();
  714. if (Receivers == nullptr) throw NullReceivers();
  715. #ifdef LEMMAUSEOMP
  716. #pragma omp parallel
  717. #endif
  718. { // OpenMP Parallel Block
  719. #ifdef LEMMAUSEOMP
  720. int tid = omp_get_thread_num();
  721. int nthreads = omp_get_num_threads();
  722. #else
  723. int tid=0;
  724. int nthreads=1;
  725. #endif
  726. auto tDipole = Dipole->Clone();
  727. std::shared_ptr<HankelTransform> Hankel;
  728. switch (HankelType) {
  729. case ANDERSON801:
  730. Hankel = FHTAnderson801::NewSP();
  731. break;
  732. case CHAVE:
  733. Hankel = GQChave::NewSP();
  734. break;
  735. case FHTKEY201:
  736. Hankel = FHTKey201::NewSP();
  737. break;
  738. case FHTKEY101:
  739. Hankel = FHTKey101::NewSP();
  740. break;
  741. case FHTKEY51:
  742. Hankel = FHTKey51::NewSP();
  743. break;
  744. case QWEKEY:
  745. Hankel = QWEKey::NewSP();
  746. break;
  747. default:
  748. std::cerr << "Hankel transform cannot be created\n";
  749. exit(EXIT_FAILURE);
  750. }
  751. if ( tDipole->GetNumberOfFrequencies() < Receivers->GetNumberOfPoints() ) {
  752. for (int ifreq=0; ifreq<tDipole->GetNumberOfFrequencies(); ++ifreq) {
  753. // Propogation constant in free space being input to Hankel
  754. Real wavef = tDipole->GetAngularFrequency(ifreq) * std::sqrt(MU0*EPSILON0);
  755. for (int irec=tid; irec<Receivers->GetNumberOfPoints(); irec+=nthreads) {
  756. SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
  757. }
  758. }
  759. } else {
  760. for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
  761. for (int ifreq=tid; ifreq<tDipole->GetNumberOfFrequencies(); ifreq+=nthreads) {
  762. // Propogation constant in free space being input to Hankel
  763. Real wavef = tDipole->GetAngularFrequency(ifreq) * std::sqrt(MU0*EPSILON0);
  764. SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
  765. }
  766. }
  767. }
  768. } // OpenMP Parallel Block
  769. }
  770. NullReceivers::NullReceivers() :
  771. runtime_error("nullptr RECEIVERS") {}
  772. NullAntenna::NullAntenna() :
  773. runtime_error("nullptr ANTENNA") {}
  774. NullInstrument::NullInstrument(LemmaObject* ptr) :
  775. runtime_error("nullptr INSTRUMENT") {
  776. std::cout << "Thrown by instance of "
  777. << ptr->GetName() << std::endl;
  778. }
  779. DipoleSourceSpecifiedForWireAntennaCalc::
  780. DipoleSourceSpecifiedForWireAntennaCalc() :
  781. runtime_error("DIPOLE SOURCE SPECIFIED FOR WIRE ANTENNA CALC"){}
  782. } // end of Lemma Namespace