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.

QWEKey.h 6.5KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202
  1. /* This file is part of Lemma, a geophysical modelling and inversion API.
  2. * More information is available at http://lemmasoftware.org
  3. */
  4. /* This Source Code Form is subject to the terms of the Mozilla Public
  5. * License, v. 2.0. If a copy of the MPL was not distributed with this
  6. * file, You can obtain one at http://mozilla.org/MPL/2.0/.
  7. */
  8. /**
  9. * @file
  10. * @date 02/12/2014 10:20:18 AM
  11. * @version $Id$
  12. * @author Trevor Irons (ti)
  13. * @email Trevor.Irons@lemmasoftware.org
  14. * @copyright Copyright (c) 2014, Trevor Irons
  15. */
  16. #ifndef QWEKEY_INC
  17. #define QWEKEY_INC
  18. #include "HankelTransform.h"
  19. #include <Eigen/Eigenvalues>
  20. #include "boost/math/special_functions.hpp"
  21. #include "boost/math/special_functions/bessel.hpp"
  22. namespace Lemma {
  23. /** breakpoint to use in division of domain, based on zeros of bessel function or
  24. regular nPi spacing.
  25. */
  26. enum sZeroType{J0, J1, NPI};
  27. /**
  28. \brief Port of Key's quadrature with extrapolation Hankel transform algorithm.
  29. \details Details of the algorithm can be found in Key2011. This code is a port
  30. of the published algorithm, which contains the following notice:
  31. %------------------------------------------------------------------%
  32. % Copyright (c) 2012 by the Society of Exploration Geophysicists. %
  33. % For more information, go to http://software.seg.org/2012/0003 . %
  34. % You must read and accept usage terms at: %
  35. % http://software.seg.org/disclaimer.txt before use. %
  36. %------------------------------------------------------------------%
  37. */
  38. class QWEKey : public HankelTransform {
  39. friend std::ostream &operator<<(std::ostream &stream, const QWEKey &ob);
  40. struct ctor_key {};
  41. public:
  42. // ==================== LIFECYCLE =======================
  43. /** Default locked constructor, use NewSP */
  44. QWEKey ( const ctor_key& );
  45. /** DeSerializing locked constructor, use DeSerialize */
  46. QWEKey ( const YAML::Node& node, const ctor_key& );
  47. /** Default destructor */
  48. ~QWEKey ();
  49. /**
  50. * Factory method for generating objects.
  51. * @return std::shared_ptr< QWEKey >
  52. */
  53. static std::shared_ptr<QWEKey> NewSP();
  54. /** YAML Serializing method
  55. */
  56. YAML::Node Serialize() const;
  57. /**
  58. * Constructs an object from a YAML::Node.
  59. */
  60. static std::shared_ptr< QWEKey > DeSerialize(const YAML::Node& node);
  61. // ==================== OPERATORS =======================
  62. void TestPrivate(const int& N);
  63. // ==================== OPERATIONS =======================
  64. Complex Zgauss(const int &ikk, const EMMODE &imode,
  65. const int &itype, const Real &rho,
  66. const Real &wavef, KernelEM1DBase* Kernel);
  67. /// Computes related kernels, if applicable, otherwise this is
  68. /// just a dummy function.
  69. void ComputeRelated(const Real& rho, std::shared_ptr<KernelEM1DBase> Kernel);
  70. void ComputeRelated(const Real& rho, std::vector< std::shared_ptr<KernelEM1DBase> > KernelVec);
  71. void ComputeRelated(const Real& rho, std::shared_ptr<KernelEM1DManager> KernelManager);
  72. // ==================== ACCESS =======================
  73. // ==================== INQUIRY =======================
  74. /** Returns the name of the underlying class, similiar to Python's type */
  75. virtual inline std::string GetName() const {
  76. return CName;
  77. }
  78. protected:
  79. // ==================== LIFECYCLE =======================
  80. /** Calculates Gauss quadrature weights of order N on the interval -1,1
  81. Algorithm from p 129 in:
  82. Trefethen, L. N., 2000, Spectral methods in MATLAB: Society for
  83. Industrial and Applied Mathematics (SIAM), volume 10 of Software,
  84. Environments, and Tools.
  85. */
  86. void GaussQuadWeights(const int& N);
  87. /** Returns the quadrature intervals and Bessel function weights used for the
  88. QWE method.
  89. */
  90. void BesselWeights( const sZeroType& sType);
  91. /** Computes an infinite integral using the partial sum of quadrature terms
  92. accelerated by sequence extrapolation using the Shanks transformation
  93. implemented with Wynn's epsilon algorithm.
  94. */
  95. void QWE(const Real& rho);
  96. /** Calls the underlying kernel functions evaluated as necessary
  97. */
  98. void getEyKernel(const int& i, const int& idx, const Real& rho);
  99. private:
  100. // ==================== DATA MEMBERS =========================
  101. /** Relative tolerance, default is 1e-6 */
  102. Real RelTol;
  103. /** Absolute tolerance, default is 1e-24 */
  104. Real AbsTol;
  105. /** Quadrature order, higher is more accurate but more expensive. Eefault is 9 */
  106. int nQuad;
  107. /** in QWE partial integrals before Shanks recurive algorithm. Defaults to 1 */
  108. int nDelay;
  109. /** Maximum number of intervals to integrate over . Defaults to 40 */
  110. int nIntervalsMax;
  111. /** Weighing of gaussian quadrature points */
  112. VectorXr GaussWeights;
  113. /** Abscissa locations of quadrature points */
  114. VectorXr GaussAbscissa;
  115. /** Breakpoints for dividing up the global integral */
  116. VectorXr xInt;
  117. /** All quadrature points between all breakpoints */
  118. VectorXr Bx;
  119. /** J0 weights */
  120. VectorXr BJ0;
  121. /** J1 weights */
  122. VectorXr BJ1;
  123. /** array of lambda arguments */
  124. VectorXr Lambda;
  125. /** array of lambda arguments */
  126. VectorXr Intervals;
  127. MatrixXcr TS;
  128. VectorXi Tn;
  129. MatrixXcr Textrap;
  130. MatrixXr TrelErr;
  131. MatrixXr TabsErr;
  132. /** Container to hold bessel arguments */
  133. Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic > Zwork;
  134. /** Container to hold bessel arguments */
  135. Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic > Zans;
  136. /** Manager for related kernels to evaluate */
  137. std::shared_ptr<KernelEM1DManager> KernelManager;
  138. /** ASCII string representation of the class name */
  139. static constexpr auto CName = "QWEKey";
  140. }; // ----- end of class QWEKey -----
  141. } // ----- end of Lemma name -----
  142. #endif // ----- #ifndef QWEKEY_INC -----