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

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