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.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200
  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 std::string GetName() const ;
  76. protected:
  77. // ==================== LIFECYCLE =======================
  78. /** Calculates Gauss quadrature weights of order N on the interval -1,1
  79. Algorithm from p 129 in:
  80. Trefethen, L. N., 2000, Spectral methods in MATLAB: Society for
  81. Industrial and Applied Mathematics (SIAM), volume 10 of Software,
  82. Environments, and Tools.
  83. */
  84. void GaussQuadWeights(const int& N);
  85. /** Returns the quadrature intervals and Bessel function weights used for the
  86. QWE method.
  87. */
  88. void BesselWeights( const sZeroType& sType);
  89. /** Computes an infinite integral using the partial sum of quadrature terms
  90. accelerated by sequence extrapolation using the Shanks transformation
  91. implemented with Wynn's epsilon algorithm.
  92. */
  93. void QWE(const Real& rho);
  94. /** Calls the underlying kernel functions evaluated as necessary
  95. */
  96. void getEyKernel(const int& i, const int& idx, const Real& rho);
  97. private:
  98. // ==================== DATA MEMBERS =========================
  99. /** Relative tolerance, default is 1e-6 */
  100. Real RelTol;
  101. /** Absolute tolerance, default is 1e-24 */
  102. Real AbsTol;
  103. /** Quadrature order, higher is more accurate but more expensive. Eefault is 9 */
  104. int nQuad;
  105. /** in QWE partial integrals before Shanks recurive algorithm. Defaults to 1 */
  106. int nDelay;
  107. /** Maximum number of intervals to integrate over . Defaults to 40 */
  108. int nIntervalsMax;
  109. /** Weighing of gaussian quadrature points */
  110. VectorXr GaussWeights;
  111. /** Abscissa locations of quadrature points */
  112. VectorXr GaussAbscissa;
  113. /** Breakpoints for dividing up the global integral */
  114. VectorXr xInt;
  115. /** All quadrature points between all breakpoints */
  116. VectorXr Bx;
  117. /** J0 weights */
  118. VectorXr BJ0;
  119. /** J1 weights */
  120. VectorXr BJ1;
  121. /** array of lambda arguments */
  122. VectorXr Lambda;
  123. /** array of lambda arguments */
  124. VectorXr Intervals;
  125. MatrixXcr TS;
  126. VectorXi Tn;
  127. MatrixXcr Textrap;
  128. MatrixXr TrelErr;
  129. MatrixXr TabsErr;
  130. /** Container to hold bessel arguments */
  131. Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic > Zwork;
  132. /** Container to hold bessel arguments */
  133. Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic > Zans;
  134. /** Manager for related kernels to evaluate */
  135. std::shared_ptr<KernelEM1DManager> KernelManager;
  136. /** ASCII string representation of the class name */
  137. static constexpr auto CName = "QWEKey";
  138. }; // ----- end of class QWEKey -----
  139. } // ----- end of Lemma name -----
  140. #endif // ----- #ifndef QWEKEY_INC -----