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

EMEarth1D.cpp 37KB

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 || IRONS )) {
  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