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.

KernelEM1DReflBase.h 7.6KB

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