Lemma is an Electromagnetics API
Você não pode selecionar mais de 25 tópicos Os tópicos devem começar com uma letra ou um número, podem incluir traços ('-') e podem ter até 35 caracteres.

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. #ifdef HAVE_BOOST_SPECIAL_FUNCTIONS
  21. #include "boost/math/special_functions.hpp"
  22. #include "boost/math/special_functions/bessel.hpp"
  23. #endif
  24. namespace Lemma {
  25. /** breakpoint to use in division of domain, based on zeros of bessel function or
  26. regular nPi spacing.
  27. */
  28. enum sZeroType{J0, J1, NPI};
  29. /**
  30. \ingroup FDEM1D
  31. \brief Port of Key's quadrature with extrapolation Hankel transform algorithm.
  32. \details Details of the algorithm can be found in Key2011. This code is a port
  33. of the published algorithm, which contains the following notice:
  34. %------------------------------------------------------------------%
  35. % Copyright (c) 2012 by the Society of Exploration Geophysicists. %
  36. % For more information, go to http://software.seg.org/2012/0003 . %
  37. % You must read and accept usage terms at: %
  38. % http://software.seg.org/disclaimer.txt before use. %
  39. %------------------------------------------------------------------%
  40. */
  41. class QWEKey : public HankelTransform {
  42. friend std::ostream &operator<<(std::ostream &stream, const QWEKey &ob);
  43. public:
  44. // ==================== LIFECYCLE =======================
  45. /** Default locked constructor, use NewSP */
  46. QWEKey ( const ctor_key& );
  47. /** DeSerializing locked constructor, use DeSerialize */
  48. QWEKey ( const YAML::Node& node, const ctor_key& );
  49. /** Default destructor */
  50. ~QWEKey ();
  51. /**
  52. * Factory method for generating objects.
  53. * @return std::shared_ptr< QWEKey >
  54. */
  55. static std::shared_ptr<QWEKey> NewSP();
  56. /** YAML Serializing method
  57. */
  58. YAML::Node Serialize() const;
  59. /**
  60. * Constructs an object from a YAML::Node.
  61. */
  62. static std::shared_ptr< QWEKey > DeSerialize(const YAML::Node& node);
  63. // ==================== OPERATORS =======================
  64. void TestPrivate(const int& N);
  65. // ==================== OPERATIONS =======================
  66. Complex Zgauss(const int &ikk, const EMMODE &imode,
  67. const int &itype, const Real &rho,
  68. const Real &wavef, KernelEM1DBase* Kernel);
  69. /// Computes related kernels, if applicable, otherwise this is
  70. /// just a dummy function.
  71. void ComputeRelated(const Real& rho, std::shared_ptr<KernelEM1DBase> Kernel);
  72. void ComputeRelated(const Real& rho, std::vector< std::shared_ptr<KernelEM1DBase> > KernelVec);
  73. void ComputeRelated(const Real& rho, std::shared_ptr<KernelEM1DManager> KernelManager);
  74. // ==================== ACCESS =======================
  75. // ==================== INQUIRY =======================
  76. /** Returns the name of the underlying class, similiar to Python's type */
  77. virtual std::string GetName() const ;
  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 -----