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.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243
  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 05/18/2012
  9. **/
  10. #ifndef KERNELEM1DREFLBASE_INC
  11. #define KERNELEM1DREFLBASE_INC
  12. #include "KernelEM1DBase.h"
  13. #include "DipoleSource.h"
  14. #include "LayeredEarthEM.h"
  15. namespace Lemma {
  16. enum DIPOLE_LOCATION { INAIR, INGROUND };
  17. // forward declaration for friend
  18. template<EMMODE Mode, int Ikernel, DIPOLE_LOCATION Isource, DIPOLE_LOCATION Irecv>
  19. class KernelEM1DSpec;
  20. // ===================================================================
  21. // Class: KernelEM1DReflBase
  22. /**
  23. \ingroup FDEM1D
  24. \brief Abstract class defining EM1DRefl class.
  25. \details Derived classes are template specialized for optimal performance.
  26. */
  27. // ===================================================================
  28. class KernelEM1DReflBase : public LemmaObject {
  29. template<EMMODE Mode, int Ikernel, DIPOLE_LOCATION Isource, DIPOLE_LOCATION Irecv>
  30. friend class KernelEM1DSpec;
  31. friend class KernelEM1DManager;
  32. //friend class DipoleSource;
  33. public:
  34. // ==================== LIFECYCLE =======================
  35. // ==================== OPERATORS =======================
  36. // ==================== OPERATIONS =======================
  37. void Initialise( std::shared_ptr<LayeredEarthEM> EarthIn) {
  38. nlay = EarthIn->GetNumberOfLayers();
  39. zh = VectorXcr::Zero(nlay);
  40. yh = VectorXcr::Zero(nlay);
  41. u = VectorXcr::Zero(nlay);
  42. cf = VectorXcr::Zero(nlay); // nlay -1 (lay 0 empty)
  43. rtu = VectorXcr::Zero(nlay); // nlay -1 Interfaces only
  44. rtd = VectorXcr::Zero(nlay); // nlay -1 Interfaces only
  45. kk = VectorXcr::Zero(nlay);
  46. Zyu = VectorXcr::Zero(nlay);
  47. Zyd = VectorXcr::Zero(nlay);
  48. Zyi = VectorXcr::Zero(nlay);
  49. th = VectorXcr::Zero(nlay);
  50. Earth = EarthIn;
  51. LayerThickness.resize(nlay);
  52. for (int ilay=0; ilay<nlay; ++ilay) {
  53. LayerThickness(ilay) = Earth->GetLayerThickness(ilay);
  54. }
  55. LayerDepth.resize(nlay);
  56. for (int ilay=0; ilay<nlay; ++ilay) {
  57. LayerDepth(ilay) = Earth->GetLayerDepth(ilay);
  58. }
  59. }
  60. void SetUpSource( DipoleSource* Dipole, const int &ifreq ) {
  61. zh(0) = Complex(0, Dipole->GetAngularFrequency(ifreq)*MU0);
  62. yh(0) = Complex(0, Dipole->GetAngularFrequency(ifreq)*EPSILON0);
  63. kk(0) = -zh(0) * yh(0);
  64. Earth->EvaluateColeColeModel(Dipole->GetAngularFrequency(ifreq));
  65. for (int ilay=1; ilay<Earth->GetNumberOfLayers(); ++ilay) {
  66. zh(ilay) = zh(0) * Earth->GetLayerSusceptibility(ilay);
  67. yh(ilay) = Earth->GetLayerConductivity(ilay) +
  68. yh(0)*Earth->GetLayerPermitivity(ilay);
  69. kk(ilay) = -zh(ilay)*yh(ilay);
  70. }
  71. tx_z = Dipole->GetLocation(2);
  72. lays = 0;
  73. Real Depth(0);
  74. for (int ilay=1; ilay<nlay; ++ilay) {
  75. if (tx_z > Depth) {
  76. lays = ilay;
  77. }
  78. Depth += LayerThickness(ilay);
  79. }
  80. }
  81. void SetUpReceiver(const Real &rz) {
  82. rx_z = rz;
  83. Real Depth(0.);
  84. layr = 0;
  85. for (int ilay=1; ilay<Earth->GetNumberOfLayers(); ++ilay) {
  86. if (rx_z > Depth) {
  87. layr = ilay;
  88. }
  89. Depth += LayerThickness(ilay);
  90. }
  91. }
  92. // ==================== ACCESS =======================
  93. Complex GetYm() {
  94. return yh(layr);
  95. }
  96. Complex GetZm() {
  97. return zh(layr);
  98. }
  99. Complex GetZs() {
  100. return zh(lays);
  101. }
  102. Complex GetKs() {
  103. return kk(lays);
  104. }
  105. // ==================== INQUIRY =======================
  106. virtual std::string GetName() const {
  107. return CName;
  108. }
  109. protected:
  110. // ==================== LIFECYCLE =======================
  111. /// Default protected constructor.
  112. KernelEM1DReflBase ( const ctor_key& key ) : LemmaObject( key )
  113. {
  114. }
  115. /// Default protected constructor.
  116. ~KernelEM1DReflBase () {
  117. }
  118. // ==================== OPERATIONS =======================
  119. /** Computes reflection coefficients. Depending on the
  120. * specialisation, this will either be TM or TE mode.
  121. */
  122. virtual void ComputeReflectionCoeffs(const Real& lambda)=0;
  123. /** Precomputes expensive calculations that are reused by insances
  124. * of KernelEM1DSpec in the calculation of Related potentials.
  125. */
  126. virtual void PreComputePotentialTerms()=0;
  127. // ==================== DATA MEMBERS =========================
  128. /// Bessel order, only 0 or 1 supported
  129. int id; // Needs to be dim nRel, or separate
  130. /// Layer the source is in
  131. int lays;
  132. /// Layer the receiver is in
  133. int layr;
  134. /// Number of Layers
  135. int nlay;
  136. /// Number of Related kernels to be computed
  137. int nRelated;
  138. /// Used in related kernel precompute calls
  139. int relIud;
  140. /// Receiver z position
  141. Real rx_z;
  142. /// Transmitter z position
  143. Real tx_z;
  144. /// bessel arg squared
  145. Real rams;
  146. /** Related part of con term */
  147. Complex relCon;
  148. Complex relenukadz;
  149. Complex relexp_pbs1;
  150. Complex relexp_pbs2;
  151. Complex rel_a;
  152. /// TM or TE mode
  153. EMMODE mode;
  154. /// Pointer to layered earth
  155. std::shared_ptr<LayeredEarthEM> Earth = nullptr;
  156. Complex uk;
  157. Complex um;
  158. VectorXcr cf; // nlay
  159. VectorXcr u; // nlay
  160. VectorXcr yh; // nlay
  161. VectorXcr zh; // nlay
  162. VectorXcr kk; // nlay
  163. VectorXcr Zyu; //(nlay); //Zyu.setZero();
  164. VectorXcr Zyd; //(nlay); //Zyd.setZero();
  165. VectorXcr Zyi; //(nlay); //Zyi.setZero();
  166. VectorXcr th; //(nlay);
  167. VectorXr LayerThickness;
  168. VectorXr LayerDepth;
  169. /// Reflection/Transmission coef. for upgoing waves in a
  170. /// layered earth model
  171. VectorXcr rtu;
  172. /// Reflection/Transmission coef. for downgoing waves in
  173. /// a layered earth model
  174. VectorXcr rtd;
  175. private:
  176. static constexpr auto CName = "KernelEM1DReflBase";
  177. }; // ----- end of class KernelEM1DReflBase -----
  178. } // ----- end of Lemma name -----
  179. #endif // ----- #ifndef KERNELEM1DREFLBASE_INC -----