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.

GQChave.h 15KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366
  1. /* This Source Code Form is subject to the terms of the Mozilla Public
  2. * License, v. 2.0. If a copy of the MPL was not distributed with this
  3. * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
  4. /**
  5. @file
  6. @author Trevor Irons
  7. @date 01/02/2010
  8. @version $Id: hankeltransformgaussianquadrature.h 199 2014-12-29 19:25:20Z tirons $
  9. **/
  10. #ifndef _HANKELTRANSFORMGAUSSIANQUADRATURE_h_INC
  11. #define _HANKELTRANSFORMGAUSSIANQUADRATURE_h_INC
  12. #include "HankelTransform.h"
  13. #include "KernelEM1DBase.h"
  14. #ifdef HAVE_BOOST_SPECIAL_FUNCTIONS
  15. #include "boost/math/special_functions.hpp"
  16. #endif
  17. namespace Lemma {
  18. // =======================================================================
  19. // Class: GQChave
  20. /// \ingroup FDEM1D
  21. /// \brief Calculates hankel transform using gaussian quadrature.
  22. /// \details Accurate but slow, this is a port of Alan Chave's public domain
  23. /// fortran code
  24. // =======================================================================
  25. class GQChave : public HankelTransform {
  26. EIGEN_MAKE_ALIGNED_OPERATOR_NEW
  27. friend std::ostream &operator<<(std::ostream &stream, const GQChave &ob);
  28. public:
  29. // ==================== LIFECYCLE ===========================
  30. /// Default locked constructor.
  31. explicit GQChave ( const ctor_key& );
  32. /** DeSerializing locked constructor, use DeSerialize */
  33. GQChave ( const YAML::Node& node, const ctor_key& );
  34. /// Default destructor
  35. virtual ~GQChave ();
  36. /**
  37. * Returns shared_ptr to new GQChave.
  38. * Location is
  39. * initialized to (0,0,0) type and polarization are
  40. * initialized to nonworking values that will throw
  41. * exceptions if used.
  42. */
  43. static std::shared_ptr<GQChave> NewSP();
  44. /** YAML Serializing method
  45. */
  46. YAML::Node Serialize() const;
  47. /**
  48. * Constructs an object from a YAML::Node.
  49. */
  50. static std::shared_ptr< GQChave > DeSerialize(const YAML::Node& node);
  51. // ==================== OPERATORS ===========================
  52. // ==================== OPERATIONS ===========================
  53. /// Performs numerical integration using Gaussian quadrature
  54. /// ikk: type of kernel depending on source and receiver couple
  55. /// imode: a switch for TE(0) and TM(1) mode
  56. /// itype: order of Bessel function
  57. /// rho is argument to integral
  58. /// wavef is the propogation constant of free space
  59. /// = omega * sqrt( EP*AMU ) amu = 4 pi e-7 ep = 8.85e-12
  60. //template <EMMODE T>
  61. Complex Zgauss(const int &ikk, const EMMODE &imode,
  62. const int &itype, const Real &rho,
  63. const Real &wavef, KernelEM1DBase* Kernel);
  64. // ==================== ACCESS ============================
  65. // ==================== INQUIRY ============================
  66. /** Returns the name of the underlying class, similiar to Python's type */
  67. virtual inline std::string GetName() const ;
  68. // ==================== DATA MEMBERS ============================
  69. protected:
  70. // ==================== OPERATIONS ============================
  71. /// Modified by Yoonho Song to branch cut, June, 1996
  72. /// Separate Gaussian quarature integral by two interval
  73. /// first: integal from 0 to wavenumber of free space
  74. /// second: integral from wavenunmber of free space to infinity
  75. /// for large arguments, it uses continued fraction also
  76. /// It is recommended to use nl = 1 to 6, nu =7
  77. /// PERFORMS AUTOMATIC CALCULATION OF BESSEL TRANSFORM TO SPECIFIED
  78. /// RELATIVE AND ABSOLUTE ERROR
  79. ///
  80. /// ARGUMENT LIST:
  81. ///
  82. /// BESR,BESI-REAL AND IMAGINARY PARTS RETURNED BY BESAUX
  83. /// iorder-ORDER OF THE BESSEL FUNCTION
  84. /// NL-LOWER LIMIT FOR GAUSS ORDER TO START COMPUTATION
  85. /// NU-UPPER LIMIT FOR GAUSS ORDER
  86. /// NU,NL=1,...7 SELECTS 3,7,15,31,63,127,AND 255 POINT GAUSS
  87. /// QUADRATURE BETWEEN THE ZERO CROSSINGS OF THE BESSEL FUNCTION
  88. /// R-ARGUMENT OF THE BESSEL FUNCTION
  89. /// RERR,AERR-RELATIVE AND ABSOLUTE ERROR FOR TERMINATION
  90. /// BESAUX TERMINATES WHEN INCREASING THE GAUSS ORDER DOES NOT
  91. /// CHANGE THE RESULT BY MORE THAN RERR OR WHEN THE ABSOLUTE ERROR
  92. /// IS LESS THAN AERR OR WHEN A GAUSS ORDER OF NU IS REACHED.
  93. /// NPCS-NUMBER OF PIECES INTO WHICH EACH PARTIAL INTEGRAND
  94. /// IS DIVIDED,
  95. /// ORDINARILY SET TO ONE. FOR VERY SMALL VALUES OF R WHERE
  96. /// THE KERNEL FUNCTION IS APPRECIABLE ONLY OVER THE FIRST FEW
  97. /// LOOPS OF THE BESSEL FUNCTION, NPCS MAY BE INCREASED TO ACHIEVE
  98. /// REASONABLE ACCURACY.
  99. /// NEW IF NEW=1, THE INTEGRANDS ARE COMPUTED AND SAVED AT EACH
  100. /// GAUSS
  101. /// ORDER. IF NEW=2, PREVIOUSLY COMPUTED INTEGRANDS ARE USED. NOTE
  102. /// THAT ORDER,R, AND NPCS MUST NOT BE CHANGED WHEN SETTING NEW=2.
  103. /// IERR-ERROR PARAMETER
  104. /// IERR=0--NORMAL RETURN
  105. /// IERR=1--RESULT NOT ACCURATE TO RERR DUE TO TOO LOW A GAUSS
  106. /// ORDER OR CONVERGENCE NOT ACHIEVED IN BESTRX
  107. //template <EMMODE T>
  108. void Besautn(Real &besr, Real &besi, const int &iorder,
  109. const int &nl, const int &nu, const Real &rho,
  110. const Real &rerr, const Real &aerr,
  111. const int &npcs, int &inew, const Real &aorb,
  112. KernelEM1DBase* Kernel);
  113. /// COMPUTES BESSEL TRANSFORM OF SPECIFIED ORDER DEFINED AS
  114. /// INTEGRAL(FUNCT(X)*J-SUB-ORDER(X*R)*DX) FROM X=0 TO INFINITY
  115. /// COMPUTATION IS ACHIEVED BY INTEGRATION BETWEEN THE ASYMPTOTIC
  116. /// ZERO CROSSINGS OF THE BESSEL FUNCTION USING GAUSS QUADRATURE.
  117. /// THE RESULTING SERIES OF PARTIAL INTEGRANDS IS SUMMED BY
  118. /// CALCULATING THE PADE APPROXIMANTS TO SPEED UP CONVERGENCE.
  119. /// ARGUMENT LIST:
  120. /// BESR,BESI REAL AND IMAGINARY PARTS RETURNED BY BESTRN
  121. /// iorder ORDER OF THE BESSEL FUNCTIONC NG NUMBER OF GAUSS
  122. /// POINTS TO USE IN THE QUADRATURE ROUTINE.
  123. /// NG=1 THROUGH 7 SELECTS 3,7,15,31,63,126,AND 255 TERMS.
  124. /// R ARGUMENT OF THE BESSEL FUNCTION
  125. /// RERR,AERR SPECIFIED RELATIVE AND ABSOLUTE ERROR FOR THE
  126. /// CALCULATION. THE INTEGRATION
  127. /// TERMINATES WHEN AN ADDITIONAL TERM DOES NOT CHANGE THE
  128. /// RESULT BY MORE THAN RERR*RESULT+AERR
  129. /// NPCS NUMBER OF PIECES INTO WHICH EACH PARTIAL I
  130. /// NTEGRAND IS DIVIDED,
  131. /// ORDINARILY SET TO ONE. FOR VERY SMALL VALUES OF RANGE
  132. /// WHERE THE KERNEL FUNCTION IS APPRECIABLE ONLY OVER THE
  133. /// FIRST FEW LOOPS OF THE BESSEL FUNCTION, NPCS MAY BE
  134. /// INCREASED TO ACHIEVE REASONABLE ACCURACY. NOTE THAT
  135. /// NPCS AFFECTS ONLY THE PADE
  136. /// SUM PORTION OF THE INTEGRATION, OVER X(NSUM) TO INFINITY.
  137. /// XSUM VECTOR OF VALUES OF THE KERNEL ARGUMENT OF FUNCT FOR WHICH
  138. /// EXPLICIT CALCULATION OF THE INTEGRAL IS DESIRED, SO THAT THE
  139. /// INTEGRAL OVER 0 TO XSUM(NSUM) IS ADDED TO THE INTEGRAL OVER
  140. /// XSUM(NSUM) TO INFINITY WITH THE PADE METHOD INVOKED ONLY FOR
  141. /// THE LATTER. THIS ALLOWS THE PADE SUMMATION METHOD TO BE
  142. /// OVERRIDDEN AND SOME TYPES OF SINGULARITIES TO BE HANDLED.
  143. /// NSUM NUMBER OF VALUES IN XSUM, MAY BE ZERO.
  144. /// NEW DETERMINES METHOD OF KERNEL CALCULATION
  145. /// NEW=0 MEANS CALCULATE BUT DO NOT SAVE INTEGRANDS
  146. /// NEW=1 MEANS CALCULATE KERNEL BY CALLING FUNCT-SAVE KERNEL
  147. /// TIMES BESSEL FUNCTION
  148. /// NEW=2 MEANS USE SAVED KERNELS TIMES BESSEL FUNCTIONS IN
  149. /// COMMON /BESINT/. NOTE THAT ORDER,R,NPCS,XSUM, AND
  150. /// NSUM MAY NOT BE CHANGED WHEN SETTING NEW=2.
  151. /// IERR ERROR PARAMETER
  152. /// 0 NORMAL RETURN-INTEGRAL CONVERGED
  153. /// 1 MEANS NO CONVERGENCE AFTER NSTOP TERMS IN THE PADE SUM
  154. ///
  155. /// SUBROUTINES REQUIRED:
  156. /// BESQUD,PADECF,CF,ZEROJ,DOT,JBESS
  157. /// A.CHAVE IGPP/UCSD
  158. /// NTERM IS MAXIMUM NUMBER OF BESSEL FUNCTION LOOPS STORED IF
  159. /// NEW.NE.0
  160. /// NSTOP IS MAXIMUM Number of Pade terms
  161. //template <EMMODE T>
  162. void Bestrn( Real &BESR, Real &BESI, const int &iorder,
  163. const int &NG, const Real &R,
  164. const Real &RERR, const Real &AERR, const int &npcs,
  165. const VectorXi& XSUM, int &NSUM, int &NEW,
  166. int &IERR, int &NCNTRL, const Real &AORB,
  167. KernelEM1DBase* Kernel);
  168. /// CALCULATES THE INTEGRAL OF F(X)*J-SUB-N(X*R) OVER THE
  169. /// INTERVAL A TO B AT A SPECIFIED GAUSS ORDER THE RESULT IS
  170. /// OBTAINED USING A SEQUENCE OF 1, 3, 7, 15, 31, 63, 127, AND 255
  171. /// POINT INTERLACING GAUSS FORMULAE SO THAT NO INTEGRAND
  172. /// EVALUATIONS ARE WASTED. THE KERNEL FUNCTIONS MAY BE
  173. /// SAVED SO THAT BESSEL TRANSFORMS OF SIMILAR KERNELS ARE COMPUTED
  174. /// WITHOUT NEW EVALUATION OF THE KERNEL. DETAILS ON THE FORMULAE
  175. /// ARE GIVEN IN 'THE OPTIMUM ADDITION OF POINTS TO QUADRATURE
  176. /// FORMULAE' BY T.N.L. PATTERSON, MATHS.COMP. 22,847-856 (1968).
  177. /// GAUSS WEIGHTS TAKEN FROM COMM. A.C.M. 16,694-699 (1973)
  178. /// ARGUMENT LIST:
  179. /// A LOWER LIMIT OF INTEGRATION
  180. /// B UPPER LIMIT OF INTEGRATION
  181. /// BESR,BESI RETURNED INTEGRAL VALUE REAL AND IMAGINARY PARTS
  182. /// NG NUMBER OF POINTS IN THE GAUSS FORMULA. NG=1,...7
  183. /// SELECTS 3,7,15,31,63,127,AND 255 POINT QUADRATURE.
  184. /// NEW SELECTS METHOD OF KERNEL EVALUATION
  185. /// NEW=0 CALCULATES KERNELS BY CALLING F - NOTHING SAVED
  186. /// NEW=1 CALCULATES KERNELS BY CALLING F AND SAVES KERNEL TIMES
  187. /// BESSEL FUNCTION IN COMMON /BESINT/
  188. /// NEW=2 USES SAVED KERNEL TIMES BESSEL FUNCTIONS IN
  189. /// COMMON /BESINT/
  190. /// iorder ORDER OF THE BESSEL FUNCTION
  191. /// R ARGUMENT OF THE BESSEL FUNCTION
  192. /// F F(X) IS THE EXTERNAL INTEGRAND SUBROUTINE
  193. /// A.CHAVE IGPP/UCSDC
  194. /// MAXIMUM NUMBER OF BESSEL FUNCTION LOOPS THAT CAN BE SAVED
  195. //template <EMMODE T>
  196. void Besqud(const Real &A, const Real &B, Real &BESR, Real &BESI,
  197. const int &NG, const int &NEW, const int &iorder,
  198. const Real &R, KernelEM1DBase* Kernel);
  199. /// COMPUTES SUM(S(I)),I=1,...N BY COMPUTATION OF PADE APPROXIMANT
  200. /// USING CONTINUED FRACTION EXPANSION. FUNCTION IS DESIGNED TO BE
  201. /// CALLED SEQUENTIALLY AS N IS INCREMENTED FROM 1 TO ITS FINAL
  202. /// VALUE. THE NTH CONTINUED FRACTION COEFFICIENT IS CALCULATED AND
  203. /// STORED AND THE NTH CONVERGENT RETURNED. IT IS UP TO THE USER TO
  204. /// STOP THE CALCULATION WHEN THE DESIRED ACCURACY IS ACHIEVED.
  205. /// ALGORITHM FROM HANGGI ET AL., Z.NATURFORSCH. 33A,402-417 (1977)
  206. /// IN THEIR NOTATION, VECTORS CFCOR,CFCOI ARE LOWER CASE D,
  207. /// VECTORS DR, DI ARE UPPER CASE D, VECTORS XR,XI ARE X, AND
  208. /// VECTORS SR,SI ARE S
  209. /// A.CHAVE IGPP/UCSD
  210. void Padecf(Real &SUMR, Real &SUMI, const int &N);
  211. /// EVALUATES A COMPLEX CONTINUED FRACTION BY RECURSIVE DIVISION
  212. /// STARTING AT THE BOTTOM, AS USED BY PADECF
  213. /// RESR,RESI ARE REAL AND IMAGINARY PARTS RETURNED
  214. /// CFCOR,CFCOI ARE REAL AND IMAGINARY VECTORS OF CONTINUED FRACTION
  215. /// COEFFICIENTS
  216. void CF( Real& RESR, Real &RESI,
  217. const VectorXr& CFCOR,
  218. const VectorXr& CFCOI,
  219. const int &N);
  220. /// COMPUTES ZERO OF BESSEL FUNCTION OF THE FIRST KIND FROM
  221. /// MCMAHON'S ASYMPTOTIC EXPANSION
  222. /// NZERO-NUMBER OF THE ZERO
  223. /// iorder-ORDER OF THE BESSEL FUNCTION (0 OR 1)
  224. Real ZeroJ(const int &ZERO, const int &IORDER);
  225. /// COMPUTES BESSEL FUNCTION OF ORDER "ORDER" AND ARGUMENT X BY
  226. /// CALLING NBS ROUTINES J0X AND J1X (REAL*8 BUT APPROXIMATELY
  227. /// REAL*4 ACCURACY).
  228. /// FOR MORE ACCURACY JBESS COULD BE CHANGED TO CALL, FOR EXAMPLE,
  229. /// THE IMSL ROUTINES MMBSJ0,MMBSJ1 << SEE C// BELOW >>
  230. Real Jbess(const Real &X, const int &IORDER);
  231. /// COMPUTES DOT PRODUCT OF TWO D.P. VECTORS WITH NONUNIT
  232. /// INCREMENTING ALLOWED. REPLACEMENT FOR BLAS SUBROUTINE SDOT.
  233. /// Currently does no checking, kind of stupid.
  234. /// The fortran version will wrap around if (inc*N) > X1.size()
  235. /// but not in a nice way.
  236. Real _dot(const int&N,
  237. const Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> &X1,
  238. const int &INC1,
  239. const Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> &X2,
  240. const int &INC2);
  241. // ==================== DATA MEMBERS ============================
  242. static const Real PI2;
  243. static const Real X01P;
  244. static const Real XMAX;
  245. static const Real XSMALL;
  246. static const Real J0_X01;
  247. static const Real J0_X02;
  248. static const Real J0_X11;
  249. static const Real J0_X12;
  250. static const Real FUDGE;
  251. static const Real FUDGEX;
  252. static const Real TWOPI1;
  253. static const Real TWOPI2;
  254. static const Real RTPI2;
  255. static const Real XMIN;
  256. static const Real J1_X01;
  257. static const Real J1_X02;
  258. static const Real J1_X11;
  259. static const Real J1_X12;
  260. /// Highest gauss order used, Was NG
  261. int HighestGaussOrder;
  262. /// Total number of partial integrals on last call, was NI
  263. int NumberPartialIntegrals;
  264. /// Total number of function calls, was NF
  265. int NumberFunctionEvals;
  266. int np;
  267. int nps;
  268. /////////////////////////////////////////////////////////////
  269. // Eigen members
  270. // Shared constant values
  271. static const VectorXr WT;
  272. static const VectorXr WA;
  273. Eigen::Matrix<int, 100, 1> Nk;
  274. //Eigen::Matrix<Real, 255, 100> karg;
  275. //Eigen::Matrix<Real, 510, 100> kern;
  276. Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> karg;
  277. Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> kern;
  278. // Was Besval COMMON block
  279. /*
  280. Eigen::Matrix<Real, 100, 1> Xr;
  281. Eigen::Matrix<Real, 100, 1> Xi;
  282. Eigen::Matrix<Real, 100, 1> Dr;
  283. Eigen::Matrix<Real, 100, 1> Di;
  284. Eigen::Matrix<Real, 100, 1> Sr;
  285. Eigen::Matrix<Real, 100, 1> Si;
  286. Eigen::Matrix<Real, 100, 1> Cfcor;
  287. Eigen::Matrix<Real, 100, 1> Cfcoi;
  288. */
  289. private:
  290. Eigen::Matrix<Real, Eigen::Dynamic, 1> Xr;
  291. Eigen::Matrix<Real, Eigen::Dynamic, 1> Xi;
  292. Eigen::Matrix<Real, Eigen::Dynamic, 1> Dr;
  293. Eigen::Matrix<Real, Eigen::Dynamic, 1> Di;
  294. Eigen::Matrix<Real, Eigen::Dynamic, 1> Sr;
  295. Eigen::Matrix<Real, Eigen::Dynamic, 1> Si;
  296. Eigen::Matrix<Real, Eigen::Dynamic, 1> Cfcor;
  297. Eigen::Matrix<Real, Eigen::Dynamic, 1> Cfcoi;
  298. /** ASCII string representation of the class name */
  299. static constexpr auto CName = "GQChave";
  300. }; // ----- end of class GQChave -----
  301. //////////////////////////////////////////////////////////////
  302. // Exception Classes
  303. /** If the lower integration limit is greater than the upper limit, throw this
  304. * error.
  305. */
  306. class LowerGaussLimitGreaterThanUpperGaussLimit :
  307. public std::runtime_error {
  308. /** Thrown when the LowerGaussLimit is greater than the upper limit.
  309. */
  310. public: LowerGaussLimitGreaterThanUpperGaussLimit();
  311. };
  312. } // ----- end of Lemma name -----
  313. #endif // ----- #ifndef _HANKELTRANSFORMGAUSSIANQUADRATURE_h_INC -----