Main Lemma Repository
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.

FHTAnderson801.h 9.8KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269
  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 06/26/2009
  9. @version $Id: hankeltransformhankel2.h 201 2015-01-03 00:07:47Z tirons $
  10. **/
  11. #ifndef __FHTANDERSON801_H
  12. #define __FHTANDERSON801_H
  13. #include "KernelEM1DBase.h"
  14. #include "KernelEM1DSpec.h"
  15. #include "CubicSplineInterpolator.h"
  16. #include "HankelTransform.h"
  17. namespace Lemma {
  18. // ==========================================================================
  19. // Class: FHTAnderson801
  20. /** \brief Computes the Hankel transform of orders 0 and 1 using lagged
  21. and related convolutions.
  22. \details A rewrite of work by Anderson who wrote a FORTRAN program
  23. that he released while working at the USGS:
  24. Anderson, W. L., 1989, A hybrid fast hankel transform algorithm for
  25. electromagnetic modeling: Geophysics, 54, 263-266.
  26. This function does not provide the Hybrid functionality however, merely the
  27. digital filter implimentation.
  28. The transform evaluates an integral of the form:
  29. \f[ \int_0^\infty K(\lambda) J_I (\lambda r) ~ d \lambda
  30. \f]
  31. Where \f$ K(\lambda) \f$ is some kernel function. The value
  32. \f$ J_I \f$ is the Bessel function of order
  33. \f$I, I \in \{0,1\} \f$
  34. The kernel function is unique for each source and is computed
  35. in the class CalculateKernel. The value \f$r\f$ is the radial
  36. distance away from the centre of the grid \f$ r=\sqrt{x^2 + y^2} \f$
  37. The Hankel transform is useful as it allows a double fourier
  38. transform to be written as a single integral:
  39. \f[ \mathop {\int \!\!\! \int}_{\!\!\!\!\!-\infty}^{\,\,\infty}
  40. F(k_x^2 + k_y^2)
  41. e^{\imath (k_x x + k_y y)} dk_x \, dk_y = 2 \pi
  42. \int_0^\infty K(\lambda) J_I (\lambda r) ~ d \lambda
  43. \f]
  44. This can only be done where there is radial symmetry. Hence
  45. its application to 1D solutions here.
  46. \note In previous versions of Lemma, this class was called HankelTransformHankel2,
  47. which more closely follows Anderson's procedural routine names, but was non-descriptive
  48. regarding where the algorithm is derived from.
  49. */
  50. // ==========================================================================
  51. class FHTAnderson801 : public HankelTransform {
  52. friend std::ostream &operator<<(std::ostream &stream, const FHTAnderson801 &ob);
  53. struct ctor_key {};
  54. public:
  55. // ==================== LIFECYCLE ==============================
  56. /**
  57. * Returns shared_ptr to new FHTAnderson801.
  58. */
  59. static std::shared_ptr<FHTAnderson801> NewSP();
  60. /**
  61. * Returns unique_ptr to new FHTAnderson801.
  62. */
  63. static std::unique_ptr<FHTAnderson801> NewUP();
  64. /// Default locked constructor
  65. FHTAnderson801( const ctor_key& );
  66. /** Locked deserializing constructor. */
  67. FHTAnderson801 ( const YAML::Node& node, const ctor_key& );
  68. /// Default destructor
  69. ~FHTAnderson801();
  70. /**
  71. * YAML Serializing method
  72. */
  73. YAML::Node Serialize() const;
  74. /**
  75. * Constructs an object from a YAML::Node.
  76. */
  77. static std::shared_ptr< FHTAnderson801 > DeSerialize(const YAML::Node& node);
  78. // ==================== OPERATORS ==============================
  79. // ==================== OPERATIONS ==============================
  80. /// Sets the number of convolutions
  81. void SetNumConv(const int &i);
  82. /// Computes the hankel transform with arguments
  83. /// @param rho [input] rho is the hankel transform argument
  84. /// @param ntol [input] ntol is
  85. /// @param tol [input] tol is
  86. void Compute(const Real &rho, const int& ntol, const Real &tol);
  87. /// Computes the related
  88. void ComputeRelated(const Real &rho, std::shared_ptr<KernelEM1DBase> Kernel);
  89. /// Computes the related
  90. void ComputeRelated(const Real &rho, std::vector< std::shared_ptr<KernelEM1DBase> > KernelVec);
  91. /// Computes the related
  92. void ComputeRelated(const Real &rho, std::shared_ptr<KernelEM1DManager> Manager);
  93. /// Computes the related and lagged convolutions
  94. void ComputeLaggedRelated(const Real &rho, const int& nlag, std::shared_ptr<KernelEM1DManager> Manager);
  95. // ==================== ACCESS ==============================
  96. /// Returns the answer
  97. Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic> GetAnswer();
  98. /// Returns the arguments for lagged convolutions
  99. VectorXr GetArg() {return Arg;};
  100. /// Returns the value of Abscissa stepping
  101. Real GetABSER( ) { return ABSER; };
  102. /// Sets the lagged kernel index so that the proper value is returned
  103. void SetLaggedArg(const Real& rho);
  104. // ==================== INQUIRY ==============================
  105. /// Calculates Hankel Transform using filtering.
  106. /// ikk: type of kernel depending on source and receiver couple
  107. /// imode: a switch for TE(0) and TM(1) mode
  108. /// itype: order of Bessel function
  109. /// rho is argument to integral
  110. /// wavef is the propogation constant of free space
  111. /// = omega * sqrt( EP*AMU ) amu = 4 pi e-7 ep = 8.85e-12
  112. Complex Zgauss(const int &ikk, const EMMODE &imode,
  113. const int &itype, const Real &rho,
  114. const Real &wavef, std::shared_ptr<KernelEM1DBase> Kernel);
  115. /** Returns the name of the underlying class, similiar to Python's type */
  116. virtual inline std::string GetName() const {
  117. return CName;
  118. }
  119. protected:
  120. private:
  121. // ==================== LIFECYCLE ==============================
  122. /** A rewrite of Anderson's "Pseudo-subroutine" computed GOTO madness. */
  123. inline void StoreRetreive(const int &idx, const int &lag,
  124. Complex &Zsum, const int &irel, Complex &C, const Real& rho0) {
  125. int look = idx+lag;
  126. int iq = look/801;
  127. int ir = look%801;
  128. int iroll = iq*800;
  129. if(this->Key[ir] <= iroll) {
  130. this->Key[ir] = iroll + ir;
  131. ++this->NumFun;
  132. Manager->ComputeReflectionCoeffs(this->Lambda, idx, rho0);
  133. for (unsigned int ir2=0; ir2<this->kernelVec.size(); ++ir2) {
  134. this->Zwork(ir, ir2) = this->kernelVec[ir2]->RelBesselArg(this->Lambda);
  135. }
  136. }
  137. C = this->Zwork(ir, irel) * this->FilterWeights(this->BesselOrder, idx);
  138. Zsum += C;
  139. return;
  140. }
  141. // ==================== OPERATIONS ==============================
  142. void DeleteSplines();
  143. // ==================== DATA MEMBERS ==============================
  144. /// The hankel transform wavenumber embedded in the integral
  145. Real Lambda;
  146. /// Number of times a kernel was evaluated
  147. int NumFun;
  148. /// Number of lagged convolutions
  149. /// must be greater or equal to 1
  150. /// It is set automatically in the @see Compute function so
  151. /// that \f$ \rho \exp\left( -.1*(\mathtt{NumConv} -1) \right) \f$
  152. /// does not underflow the exponent range
  153. int NumConv;
  154. /// Number of related kernels
  155. int NumRel;
  156. /** Bessel transform order to use */
  157. int BesselOrder;
  158. /** Lag argument */
  159. int iLag;
  160. /* Should results be cached? Useful for repeated calculations of few receiver points */
  161. // turned out to have only marginal benefits in best case, and awful consequences in many
  162. //bool cacheResults;
  163. /** Related Kernel Manager */
  164. std::shared_ptr<KernelEM1DManager> Manager;
  165. /// Used as base for filter abscissa generation
  166. static const Real ABSCISSA;
  167. /// Also used in abscissa generation \f$ ABSE = \exp{.1} \f$
  168. static const Real ABSE;
  169. /// Also used in abscissa generation \f$ ABSER = 1 / \exp{.1} \f$
  170. static const Real ABSER;
  171. /// Counter for calculated
  172. int icount;
  173. /// Kernel Calculator
  174. std::vector < std::shared_ptr<KernelEM1DBase> > kernelVec;
  175. /// Spines for lagged convolutions (real part)
  176. std::vector <std::shared_ptr<CubicSplineInterpolator> > splineVecReal;
  177. /// Spines for lagged convolutions (imaginary part)
  178. std::vector < std::shared_ptr<CubicSplineInterpolator> > splineVecImag;
  179. /// Key used internally
  180. Eigen::Matrix<int, 801, 1> Key;
  181. //int Key[801];
  182. //Eigen::Matrix<int, Eigen::Dynamic, 1> Key;
  183. /// Filter weight coefficients. Set for either \f$J_0\f$ or \f$J_1\f$
  184. /// internally by protected function SetFilterWeights.
  185. /// Fixed sized will yield better performance. (not necessarily)
  186. //Eigen::Matrix<Real, 801, 1> FilterWeights;
  187. static const Eigen::Matrix<Real, 2, 801> FilterWeights;
  188. //static const Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> FilterWeights;
  189. /// Zwork from Anderson
  190. Eigen::Matrix<Complex, 801, Eigen::Dynamic> Zwork;
  191. //Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic> Zwork;
  192. /// Holds answer, dimensions are NumConv, and NumberRelated.
  193. Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic> Zans;
  194. /// Holds the arguments for lagged convolutions
  195. VectorXr Arg;
  196. /** ASCII string representation of the class name */
  197. static constexpr auto CName = "FHTAnderson801";
  198. }; // ----- end of class FHTAnderson801 -----
  199. }
  200. #endif // __FHTAnderson801_h