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

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