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 38KB

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