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.

GQChave.h 14KB

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