Lemma is an Electromagnetics API
Du kannst nicht mehr als 25 Themen auswählen Themen müssen mit entweder einem Buchstaben oder einer Ziffer beginnen. Sie können Bindestriche („-“) enthalten und bis zu 35 Zeichen lang sein.

GQChave.cpp 42KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149
  1. /* This file is part of Lemma, a geophysical modelling and inversion API */
  2. /* This Source Code Form is subject to the terms of the Mozilla Public
  3. * License, v. 2.0. If a copy of the MPL was not distributed with this
  4. * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
  5. /**
  6. @file
  7. @author Trevor Irons
  8. @date 01/02/2010
  9. @version $Id: hankeltransformgaussianquadrature.cpp 200 2014-12-29 21:11:55Z tirons $
  10. **/
  11. // Description: Port of Alan Chave's gaussian quadrature algorithm, which
  12. // is public domain, code listed in published article:
  13. // Chave, A. D., 1983, Numerical integration of related Hankel transforms by
  14. // quadrature and continued fraction expansion: Geophysics, 48
  15. // 1671--1686 doi: 10.1190/1.1441448
  16. #include "GQChave.h"
  17. namespace Lemma{
  18. std::ostream &operator<<(std::ostream &stream, const GQChave &ob) {
  19. stream << ob.Serialize() << "\n";
  20. return stream;
  21. }
  22. // Initialise static members
  23. const VectorXr GQChave::WT =
  24. (VectorXr(254) << // (WT(I),I=1,20)
  25. 0.55555555555555555556e+00,0.88888888888888888889e+00,
  26. 0.26848808986833344073e+00,0.10465622602646726519e+00,
  27. 0.40139741477596222291e+00,0.45091653865847414235e+00,
  28. 0.13441525524378422036e+00,0.51603282997079739697e-01,
  29. 0.20062852937698902103e+00,0.17001719629940260339e-01,
  30. 0.92927195315124537686e-01,0.17151190913639138079e+00,
  31. 0.21915685840158749640e+00,0.22551049979820668739e+00,
  32. 0.67207754295990703540e-01,0.25807598096176653565e-01,
  33. 0.10031427861179557877e+00,0.84345657393211062463e-02,
  34. 0.46462893261757986541e-01,0.85755920049990351154e-01,
  35. // (WT(I),I=21,40)
  36. 0.10957842105592463824e+00,0.25447807915618744154e-02,
  37. 0.16446049854387810934e-01,0.35957103307129322097e-01,
  38. 0.56979509494123357412e-01,0.76879620499003531043e-01,
  39. 0.93627109981264473617e-01,0.10566989358023480974e+00,
  40. 0.11195687302095345688e+00,0.11275525672076869161e+00,
  41. 0.33603877148207730542e-01,0.12903800100351265626e-01,
  42. 0.50157139305899537414e-01,0.42176304415588548391e-02,
  43. 0.23231446639910269443e-01,0.42877960025007734493e-01,
  44. 0.54789210527962865032e-01,0.12651565562300680114e-02,
  45. 0.82230079572359296693e-02,0.17978551568128270333e-01,
  46. // (WT(I),I=41,60)
  47. 0.28489754745833548613e-01,0.38439810249455532039e-01,
  48. 0.46813554990628012403e-01,0.52834946790116519862e-01,
  49. 0.55978436510476319408e-01,0.36322148184553065969e-03,
  50. 0.25790497946856882724e-02,0.61155068221172463397e-02,
  51. 0.10498246909621321898e-01,0.15406750466559497802e-01,
  52. 0.20594233915912711149e-01,0.25869679327214746911e-01,
  53. 0.31073551111687964880e-01,0.36064432780782572640e-01,
  54. 0.40715510116944318934e-01,0.44914531653632197414e-01,
  55. 0.48564330406673198716e-01,0.51583253952048458777e-01,
  56. 0.53905499335266063927e-01,0.55481404356559363988e-01,
  57. // (WT(I),I=61,80)
  58. 0.56277699831254301273e-01,0.56377628360384717388e-01,
  59. 0.16801938574103865271e-01,0.64519000501757369228e-02,
  60. 0.25078569652949768707e-01,0.21088152457266328793e-02,
  61. 0.11615723319955134727e-01,0.21438980012503867246e-01,
  62. 0.27394605263981432516e-01,0.63260731936263354422e-03,
  63. 0.41115039786546930472e-02,0.89892757840641357233e-02,
  64. 0.14244877372916774306e-01,0.19219905124727766019e-01,
  65. 0.23406777495314006201e-01,0.26417473395058259931e-01,
  66. 0.27989218255238159704e-01,0.18073956444538835782e-03,
  67. 0.12895240826104173921e-02,0.30577534101755311361e-02,
  68. // (WT(I),I=81,100)
  69. 0.52491234548088591251e-02,0.77033752332797418482e-02,
  70. 0.10297116957956355524e-01,0.12934839663607373455e-01,
  71. 0.15536775555843982440e-01,0.18032216390391286320e-01,
  72. 0.20357755058472159467e-01,0.22457265826816098707e-01,
  73. 0.24282165203336599358e-01,0.25791626976024229388e-01,
  74. 0.26952749667633031963e-01,0.27740702178279681994e-01,
  75. 0.28138849915627150636e-01,0.50536095207862517625e-04,
  76. 0.37774664632698466027e-03,0.93836984854238150079e-03,
  77. 0.16811428654214699063e-02,0.25687649437940203731e-02,
  78. 0.35728927835172996494e-02,0.46710503721143217474e-02,
  79. // (WT(I),I=101,120)
  80. 0.58434498758356395076e-02,0.70724899954335554680e-02,
  81. 0.83428387539681577056e-02,0.96411777297025366953e-02,
  82. 0.10955733387837901648e-01,0.12275830560082770087e-01,
  83. 0.13591571009765546790e-01,0.14893641664815182035e-01,
  84. 0.16173218729577719942e-01,0.17421930159464173747e-01,
  85. 0.18631848256138790186e-01,0.19795495048097499488e-01,
  86. 0.20905851445812023852e-01,0.21956366305317824939e-01,
  87. 0.22940964229387748761e-01,0.23854052106038540080e-01,
  88. 0.24690524744487676909e-01,0.25445769965464765813e-01,
  89. 0.26115673376706097680e-01,0.26696622927450359906e-01,
  90. // (WT(I),I=121,140)
  91. 0.27185513229624791819e-01,0.27579749566481873035e-01,
  92. 0.27877251476613701609e-01,0.28076455793817246607e-01,
  93. 0.28176319033016602131e-01,0.28188814180192358694e-01,
  94. 0.84009692870519326354e-02,0.32259500250878684614e-02,
  95. 0.12539284826474884353e-01,0.10544076228633167722e-02,
  96. 0.58078616599775673635e-02,0.10719490006251933623e-01,
  97. 0.13697302631990716258e-01,0.31630366082226447689e-03,
  98. 0.20557519893273465236e-02,0.44946378920320678616e-02,
  99. 0.71224386864583871532e-02,0.96099525623638830097e-02,
  100. 0.11703388747657003101e-01,0.13208736697529129966e-01,
  101. // (WT(I),I=141,160)
  102. 0.13994609127619079852e-01,0.90372734658751149261e-04,
  103. 0.64476204130572477933e-03,0.15288767050877655684e-02,
  104. 0.26245617274044295626e-02,0.38516876166398709241e-02,
  105. 0.51485584789781777618e-02,0.64674198318036867274e-02,
  106. 0.77683877779219912200e-02,0.90161081951956431600e-02,
  107. 0.10178877529236079733e-01,0.11228632913408049354e-01,
  108. 0.12141082601668299679e-01,0.12895813488012114694e-01,
  109. 0.13476374833816515982e-01,0.13870351089139840997e-01,
  110. 0.14069424957813575318e-01,0.25157870384280661489e-04,
  111. 0.18887326450650491366e-03,0.46918492424785040975e-03,
  112. // (WT(I),I=161,180)
  113. 0.84057143271072246365e-03,0.12843824718970101768e-02,
  114. 0.17864463917586498247e-02,0.23355251860571608737e-02,
  115. 0.29217249379178197538e-02,0.35362449977167777340e-02,
  116. 0.41714193769840788528e-02,0.48205888648512683476e-02,
  117. 0.54778666939189508240e-02,0.61379152800413850435e-02,
  118. 0.67957855048827733948e-02,0.74468208324075910174e-02,
  119. 0.80866093647888599710e-02,0.87109650797320868736e-02,
  120. 0.93159241280693950932e-02,0.98977475240487497440e-02,
  121. 0.10452925722906011926e-01,0.10978183152658912470e-01,
  122. 0.11470482114693874380e-01,0.11927026053019270040e-01,
  123. // (WT(I),I=181,200)
  124. 0.12345262372243838455e-01,0.12722884982732382906e-01,
  125. 0.13057836688353048840e-01,0.13348311463725179953e-01,
  126. 0.13592756614812395910e-01,0.13789874783240936517e-01,
  127. 0.13938625738306850804e-01,0.14038227896908623303e-01,
  128. 0.14088159516508301065e-01,0.69379364324108267170e-05,
  129. 0.53275293669780613125e-04,0.13575491094922871973e-03,
  130. 0.24921240048299729402e-03,0.38974528447328229322e-03,
  131. 0.55429531493037471492e-03,0.74028280424450333046e-03,
  132. 0.94536151685852538246e-03,0.11674841174299594077e-02,
  133. 0.14049079956551446427e-02,0.16561127281544526052e-02,
  134. // (WT(I),I=201,220)
  135. 0.19197129710138724125e-02,0.21944069253638388388e-02,
  136. 0.24789582266575679307e-02,0.27721957645934509940e-02,
  137. 0.30730184347025783234e-02,0.33803979910869203823e-02,
  138. 0.36933779170256508183e-02,0.40110687240750233989e-02,
  139. 0.43326409680929828545e-02,0.46573172997568547773e-02,
  140. 0.49843645647655386012e-02,0.53130866051870565663e-02,
  141. 0.56428181013844441585e-02,0.59729195655081658049e-02,
  142. 0.63027734490857587172e-02,0.66317812429018878941e-02,
  143. 0.69593614093904229394e-02,0.72849479805538070639e-02,
  144. 0.76079896657190565832e-02,0.79279493342948491103e-02,
  145. // (WT(I),I=221,240)
  146. 0.82443037630328680306e-02,0.85565435613076896192e-02,
  147. 0.88641732094824942641e-02,0.91667111635607884067e-02,
  148. 0.94636899938300652943e-02,0.97546565363174114611e-02,
  149. 0.10039172044056840798e-01,0.10316812330947621682e-01,
  150. 0.10587167904885197931e-01,0.10849844089337314099e-01,
  151. 0.11104461134006926537e-01,0.11350654315980596602e-01,
  152. 0.11588074033043952568e-01,0.11816385890830235763e-01,
  153. 0.12035270785279562630e-01,0.12244424981611985899e-01,
  154. 0.12443560190714035263e-01,0.12632403643542078765e-01,
  155. 0.12810698163877361967e-01,0.12978202239537399286e-01,
  156. // (WT(I),I=241,254)
  157. 0.13134690091960152836e-01,0.13279951743930530650e-01,
  158. 0.13413793085110098513e-01,0.13536035934956213614e-01,
  159. 0.13646518102571291428e-01,0.13745093443001896632e-01,
  160. 0.13831631909506428676e-01,0.13906019601325461264e-01,
  161. 0.13968158806516938516e-01,0.14017968039456608810e-01,
  162. 0.14055382072649964277e-01,0.14080351962553661325e-01,
  163. 0.14092845069160408355e-01,0.14094407090096179347e-01).finished();
  164. const VectorXr GQChave::WA =
  165. (VectorXr(127) << // (WT(I),I=1,20)
  166. // (WA(I),I=1,20)
  167. 0.77459666924148337704e+00,0.96049126870802028342e+00,
  168. 0.43424374934680255800e+00,0.99383196321275502221e+00,
  169. 0.88845923287225699889e+00,0.62110294673722640294e+00,
  170. 0.22338668642896688163e+00,0.99909812496766759766e+00,
  171. 0.98153114955374010687e+00,0.92965485742974005667e+00,
  172. 0.83672593816886873550e+00,0.70249620649152707861e+00,
  173. 0.53131974364437562397e+00,0.33113539325797683309e+00,
  174. 0.11248894313318662575e+00,0.99987288812035761194e+00,
  175. 0.99720625937222195908e+00,0.98868475754742947994e+00,
  176. 0.97218287474858179658e+00,0.94634285837340290515e+00,
  177. // (WA(I),I=21,40)
  178. 0.91037115695700429250e+00,0.86390793819369047715e+00,
  179. 0.80694053195021761186e+00,0.73975604435269475868e+00,
  180. 0.66290966002478059546e+00,0.57719571005204581484e+00,
  181. 0.48361802694584102756e+00,0.38335932419873034692e+00,
  182. 0.27774982202182431507e+00,0.16823525155220746498e+00,
  183. 0.56344313046592789972e-01,0.99998243035489159858e+00,
  184. 0.99959879967191068325e+00,0.99831663531840739253e+00,
  185. 0.99572410469840718851e+00,0.99149572117810613240e+00,
  186. 0.98537149959852037111e+00,0.97714151463970571416e+00,
  187. 0.96663785155841656709e+00,0.95373000642576113641e+00,
  188. // (WA(I),I=41,60)
  189. 0.93832039777959288365e+00,0.92034002547001242073e+00,
  190. 0.89974489977694003664e+00,0.87651341448470526974e+00,
  191. 0.85064449476835027976e+00,0.82215625436498040737e+00,
  192. 0.79108493379984836143e+00,0.75748396638051363793e+00,
  193. 0.72142308537009891548e+00,0.68298743109107922809e+00,
  194. 0.64227664250975951377e+00,0.59940393024224289297e+00,
  195. 0.55449513263193254887e+00,0.50768775753371660215e+00,
  196. 0.45913001198983233287e+00,0.40897982122988867241e+00,
  197. 0.35740383783153215238e+00,0.30457644155671404334e+00,
  198. 0.25067873030348317661e+00,0.19589750271110015392e+00,
  199. // (WA(I),I=61,80)
  200. 0.14042423315256017459e+00,0.84454040083710883710e-01,
  201. 0.28184648949745694339e-01,0.99999759637974846462e+00,
  202. 0.99994399620705437576e+00,0.99976049092443204733e+00,
  203. 0.99938033802502358193e+00,0.99874561446809511470e+00,
  204. 0.99780535449595727456e+00,0.99651414591489027385e+00,
  205. 0.99483150280062100052e+00,0.99272134428278861533e+00,
  206. 0.99015137040077015918e+00,0.98709252795403406719e+00,
  207. 0.98351865757863272876e+00,0.97940628167086268381e+00,
  208. 0.97473445975240266776e+00,0.96948465950245923177e+00,
  209. 0.96364062156981213252e+00,0.95718821610986096274e+00,
  210. // (WA(I),I=81,100)
  211. 0.95011529752129487656e+00,0.94241156519108305981e+00,
  212. 0.93406843615772578800e+00,0.92507893290707565236e+00,
  213. 0.91543758715576504064e+00,0.90514035881326159519e+00,
  214. 0.89418456833555902286e+00,0.88256884024734190684e+00,
  215. 0.87029305554811390585e+00,0.85735831088623215653e+00,
  216. 0.84376688267270860104e+00,0.82952219463740140018e+00,
  217. 0.81462878765513741344e+00,0.79909229096084140180e+00,
  218. 0.78291939411828301639e+00,0.76611781930376009072e+00,
  219. 0.74869629361693660282e+00,0.73066452124218126133e+00,
  220. 0.71203315536225203459e+00,0.69281376977911470289e+00,
  221. // (WA(I),I=101,120)
  222. 0.67301883023041847920e+00,0.65266166541001749610e+00,
  223. 0.63175643771119423041e+00,0.61031811371518640016e+00,
  224. 0.58836243444766254143e+00,0.56590588542365442262e+00,
  225. 0.54296566649831149049e+00,0.51955966153745702199e+00,
  226. 0.49570640791876146017e+00,0.47142506587165887693e+00,
  227. 0.44673538766202847374e+00,0.42165768662616330006e+00,
  228. 0.39621280605761593918e+00,0.37042208795007823014e+00,
  229. 0.34430734159943802278e+00,0.31789081206847668318e+00,
  230. 0.29119514851824668196e+00,0.26424337241092676194e+00,
  231. 0.23705884558982972721e+00,0.20966523824318119477e+00,
  232. // (WA(I),I=121,127)
  233. 0.18208649675925219825e+00,0.15434681148137810869e+00,
  234. 0.12647058437230196685e+00,0.98482396598119202090e-01,
  235. 0.70406976042855179063e-01,0.42269164765363603212e-01,
  236. 0.14093886410782462614e-01).finished();
  237. /*
  238. const Real PI2 = 0.6366197723675813;
  239. const Real X01P = 0.4809651115391545e01;
  240. const Real XMAX = std::numeric_limits<Real>::max();
  241. const Real XSMALL = 0.9094947017729281e-12;
  242. const Real J0_X01 = 0.2404825557695772e01;
  243. const Real J0_X02 = 0.1043754397719454e-15;
  244. const Real J0_X11 = 0.5520078110286310e01;
  245. const Real J0_X12 = 0.8088597146146419e-16;
  246. const Real FUDGE = 6.071532166000000e-18;
  247. const Real FUDGEX = 1.734723476000000e-18;
  248. const Real TWOPI1 = 0.6283185005187988e01;
  249. const Real TWOPI2 = 0.3019915981956752e-06;
  250. const Real RTPI2 = 0.7978845608028652e0;
  251. const Real XMIN = std::numeric_limits<Real>::min();
  252. const Real J1_X01 = 0.3831705970207512e1;
  253. const Real J1_X02 = -0.5967810507509414e-15;
  254. const Real J1_X11 = 0.7015586669815619e1;
  255. const Real J1_X12 = -0.5382308663841630e-15;
  256. */
  257. GQChave::GQChave( const ctor_key& key ) : HankelTransform( key ) {
  258. karg.resize(255, 100);
  259. kern.resize(510, 100);
  260. }
  261. GQChave::GQChave( const YAML::Node& node, const ctor_key& key ) : HankelTransform(node, key) {
  262. }
  263. /////////////////////////////////////////////////////////////
  264. GQChave::~GQChave() {
  265. }
  266. /////////////////////////////////////////////////////////////
  267. std::shared_ptr<GQChave> GQChave::NewSP() {
  268. return std::make_shared<GQChave>( ctor_key() );
  269. }
  270. //--------------------------------------------------------------------------------------
  271. // Class: GQChave
  272. // Method: DeSerialize
  273. // Description: Factory method, converts YAML node into object
  274. //--------------------------------------------------------------------------------------
  275. std::shared_ptr<GQChave> GQChave::DeSerialize( const YAML::Node& node ) {
  276. if (node.Tag() != "GQChave") {
  277. throw DeSerializeTypeMismatch( "GQChave", node.Tag());
  278. }
  279. return std::make_shared<GQChave> ( node, ctor_key() );
  280. }
  281. //--------------------------------------------------------------------------------------
  282. // Class: GQChave
  283. // Method: Serialize
  284. // Description: Converts object into Serialized version
  285. //--------------------------------------------------------------------------------------
  286. YAML::Node GQChave::Serialize() const {
  287. YAML::Node node = HankelTransform::Serialize();
  288. node.SetTag( GetName() );
  289. //node["LayerConductivity"] = LayerConductivity;
  290. return node;
  291. }
  292. /////////////////////////////////////////////////////////////
  293. //--------------------------------------------------------------------------------------
  294. // Class: GQChave
  295. // Method: GetName
  296. // Description: Class identifier
  297. //--------------------------------------------------------------------------------------
  298. inline std::string GQChave::GetName ( ) const {
  299. return CName;
  300. } // ----- end of method GQChave::GetName -----
  301. Complex GQChave::
  302. Zgauss(const int &ikk, const EMMODE &mode,
  303. const int &itype, const Real &rho, const Real &wavef,
  304. KernelEM1DBase* Kernel){
  305. // TI, TODO, change calls to Zgauss to reflect this, go and fix so we
  306. // dont subract 1 from this everywhere
  307. //Kernel->SetIk(ikk+1);
  308. //Kernel->SetMode(mode);
  309. //ik = ikk+1;
  310. //mode = imode;
  311. Real Besr(0);
  312. Real Besi(0);
  313. // Parameters
  314. int nl(1); // Lower limit for gauss order to start comp
  315. int nu(7); // upper limit for gauss order
  316. #ifdef LEMMA_SINGLE_PRECISION
  317. Real rerr = 1e-5; // Error, for double Kihand set to .1e-10, .1e-11
  318. Real aerr = 1e-6;
  319. #else // ----- not LEMMA_SINGLE_PRECISION -----
  320. Real rerr = 1e-11; // Error, for double Kihand set to .1e-10, .1e-11
  321. Real aerr = 1e-12;
  322. #endif // ----- not LEMMA_SINGLE_PRECISION -----
  323. int npcs(1);
  324. int inew(0);
  325. //const int NTERM = 100;
  326. //BESINT.karg.resize(255, NTERM);
  327. //BESINT.kern.resize(510, NTERM);
  328. //this->karg.setZero();
  329. //this->kern.setZero();
  330. Besautn(Besr, Besi, itype, nl, nu, rho, rerr, aerr, npcs, inew, wavef, Kernel);
  331. return Complex(Besr, Besi);
  332. }
  333. //////////////////////////////////////////////////////////////////
  334. void GQChave::
  335. Besautn(Real &besr, Real &besi,
  336. const int &besselOrder,
  337. const int &lowerGaussLimit,
  338. const int &upperGaussLimit,
  339. const Real &rho,
  340. const Real &relativeError,
  341. const Real &absError,
  342. const int& numPieces,
  343. int &inew,
  344. const Real &aorb,
  345. KernelEM1DBase* Kernel) {
  346. HighestGaussOrder = 0;
  347. NumberPartialIntegrals = 0;
  348. NumberFunctionEvals = 0;
  349. inew = 0;
  350. if (lowerGaussLimit > upperGaussLimit) {
  351. besr = 0;
  352. besi = 0;
  353. throw LowerGaussLimitGreaterThanUpperGaussLimit();
  354. }
  355. int ncntrl = 0;
  356. int nw = std::max(inew, 1);
  357. // temps
  358. Real besr_1(0);
  359. Real besi_1(0);
  360. int ierr(0);
  361. int ierr1(0);
  362. int ierr2(0);
  363. VectorXi xsum(1);
  364. //xsum.setZero(); // TODO xsum doesn't do a god damn thing
  365. int nsum(0);
  366. // Check for Rtud
  367. Bestrn(besr_1, besi_1, besselOrder, lowerGaussLimit, rho,
  368. .1*relativeError, .1*absError,
  369. numPieces, xsum, nsum, nw, ierr, ncntrl, aorb, Kernel);
  370. if (ierr != 0 && lowerGaussLimit == 7) {
  371. HighestGaussOrder = lowerGaussLimit;
  372. return;
  373. } else {
  374. Real oldr = besr_1;
  375. Real oldi = besi_1;
  376. for (int n=lowerGaussLimit+1; n<=upperGaussLimit; ++n) {
  377. int nw2 = 2;
  378. Bestrn(besr_1, besi_1, besselOrder, n, rho, .1*relativeError,
  379. .1*absError, numPieces, xsum, nsum, nw2, ierr,
  380. ncntrl, aorb, Kernel);
  381. if (ierr != 0 && n==7) {
  382. besr_1 = oldr;
  383. besi_1 = oldi;
  384. HighestGaussOrder = n;
  385. std::cerr << "CONVERGENCE FAILED AT SMALL ARGUMNENT!\n";
  386. ierr1 = ierr + 10;
  387. break;
  388. } else if (std::abs(besr_1-oldr) <=
  389. relativeError*std::abs(besr_1)+absError &&
  390. std::abs(besi_1-oldi) <=
  391. relativeError*std::abs(besi_1)+absError) {
  392. HighestGaussOrder = n;
  393. break;
  394. } else {
  395. oldr = besr_1;
  396. oldi = besi_1;
  397. }
  398. }
  399. }
  400. inew = 0;
  401. ncntrl = 1;
  402. nw=std::max(inew, 1);
  403. Real besr_2, besi_2;
  404. //karg.setZero();
  405. //kern.setZero();
  406. HighestGaussOrder = 0;
  407. NumberPartialIntegrals = 0;
  408. NumberFunctionEvals = 0;
  409. Bestrn(besr_2, besi_2, besselOrder, lowerGaussLimit, rho,
  410. .1*relativeError, .1*absError,
  411. numPieces, xsum, nsum, nw, ierr, ncntrl, aorb, Kernel);
  412. if (ierr != 0 && lowerGaussLimit == 7) {
  413. HighestGaussOrder = lowerGaussLimit;
  414. return;
  415. } else {
  416. Real oldr = besr_2;
  417. Real oldi = besi_2;
  418. for (int n=lowerGaussLimit+1; n<=upperGaussLimit; ++n) {
  419. int nw2 = 2;
  420. Bestrn(besr_2, besi_2, besselOrder, n, rho, .1*relativeError,
  421. .1*absError, numPieces, xsum, nsum, nw2, ierr,
  422. ncntrl, aorb, Kernel);
  423. if (ierr != 0 && n==7) {
  424. besr_2 = oldr;
  425. besi_2 = oldi;
  426. HighestGaussOrder = n;
  427. std::cerr << "CONVERGENCE FAILED AT SMALL ARGUMNENT!\n";
  428. ierr2 = ierr + 20;
  429. break;
  430. } else if (std::abs(besr_2-oldr) <=
  431. relativeError*std::abs(besr_2)+absError &&
  432. std::abs(besi_2-oldi) <=
  433. relativeError*std::abs(besi_2)+absError) {
  434. HighestGaussOrder = n;
  435. break;
  436. } else {
  437. oldr = besr_2;
  438. oldi = besi_2;
  439. }
  440. }
  441. }
  442. besr = besr_1 + besr_2;
  443. besi = besi_1 + besi_2;
  444. ierr = ierr1 + ierr2;
  445. return;
  446. }
  447. /////////////////////////////////////////////////////////////
  448. void GQChave::
  449. Bestrn(Real &BESR, Real &BESI, const int &IORDER,
  450. const int &NG, const Real &rho,
  451. const Real &RERR, const Real &AERR, const int &NPCS,
  452. VectorXi &XSUM, int &NSUM, int &NEW,
  453. int &IERR, int &NCNTRL, const Real &AORB, KernelEM1DBase* Kernel) {
  454. Xr.setZero();
  455. Xi.setZero();
  456. Dr.setZero();
  457. Di.setZero();
  458. Sr.setZero();
  459. Si.setZero();
  460. Cfcor.setZero();
  461. Cfcoi.setZero();
  462. Dr.setZero();
  463. Di.setZero();
  464. Dr(0) = (Real)(-1);
  465. const int NTERM = 100;
  466. const int NSTOP = 100;
  467. int NPO;
  468. //std::cout << "Bestrn NEW " << NEW << std::endl;
  469. if (NEW == 2) {
  470. NPO = nps;
  471. } else {
  472. Nk.setZero();
  473. nps = 0;
  474. NPO=NTERM;
  475. }
  476. // Trivial?
  477. if (IORDER != 0 && rho == 0) {
  478. BESR = 0;
  479. BESI = 0;
  480. IERR = 0;
  481. return;
  482. }
  483. NumberPartialIntegrals=0;
  484. int NW = NEW;
  485. np = 0; // TI, zero based indexing
  486. int NPB = 1; // 0?
  487. int L = 0; // TODO, should be 0?
  488. Real B = 0.;
  489. Real A = 0.;
  490. Real SUMR = 0.;
  491. Real SUMI = 0.;
  492. Real XSUMR = 0.;
  493. Real XSUMI = 0.;
  494. Real TERMR(0), TERMI(0);
  495. // COMPUTE BESSEL TRANSFORM EXPLICITLY ON (0,XSUM(NSUM))
  496. if (NSUM > 0) {
  497. std::cerr << "NSUM GREATER THAN ZERO UNTESTED" << std::endl;
  498. Real LASTR=0.0;
  499. Real LASTI=0.0;
  500. for (int N=1; N<=NSUM; ++N) {
  501. if (NW == 2 && np > NPO) NW=1;
  502. if (np > NTERM) NW=0;
  503. A=B;
  504. B=XSUM(N);
  505. Besqud(A, B, TERMR, TERMI, NG, NW, IORDER, rho, Kernel);
  506. XSUMR += TERMR;
  507. XSUMI += TERMI;
  508. if ( (std::abs(XSUMR-LASTR) <= RERR*std::abs(XSUMR)+AERR) &&
  509. (std::abs(XSUMI-LASTI) <= RERR*std::abs(XSUMI)+AERR)) {
  510. BESR=XSUMR;
  511. BESI=XSUMI;
  512. IERR=0;
  513. nps = std::max(np, nps);
  514. return;
  515. } else {
  516. ++np;
  517. LASTR=XSUMR;
  518. LASTI=XSUMI;
  519. }
  520. }
  521. while (ZeroJ(NPB,IORDER) > XSUM(NSUM*rho)) {
  522. ++NPB;
  523. }
  524. }
  525. // ENTRY POINT FOR PADE SUMMATION OF PARTIAL INTEGRANDS
  526. Real LASTR=0.e0;
  527. Real LASTI=0.e0;
  528. if (NCNTRL == 0) {
  529. A = 0.;
  530. B = ZeroJ(NPB,IORDER) / rho;
  531. if (B > AORB) {
  532. B = AORB;
  533. }
  534. } else {
  535. A = AORB;
  536. B = ZeroJ(NPB,IORDER)/rho;
  537. while (B <= A) {
  538. ++NPB;
  539. B = ZeroJ(NPB,IORDER)/rho;
  540. }
  541. }
  542. // CALCULATE TERMS AND SUM WITH PADECF, QUITTING WHEN CONVERGENCE IS
  543. // OBTAINED
  544. if (NCNTRL != 0) {
  545. for (int N=1; N<=NSTOP; ++N) {
  546. if (NPCS == 1) {
  547. if ((NW==2) && (np > NPO)) {
  548. NW=1;
  549. }
  550. if (np > NTERM) {
  551. NW=0;
  552. }
  553. Besqud(A,B,TERMR,TERMI,NG,NW,IORDER,rho,Kernel);
  554. ++np;
  555. } else {
  556. std::cout << "In the else conditional\n";
  557. TERMR=0.;
  558. TERMI=0.;
  559. Real XINC=(B-A)/NPCS;
  560. Real AA=A;
  561. Real BB=A+XINC;
  562. for (int I=1; I<=NPCS; ++I) {
  563. if ((NW == 2) && (np > NPO)) NW=1;
  564. if (np > NTERM) NW=0;
  565. Real TR, TI;
  566. Besqud(AA, BB, TR, TI, NG, NW, IORDER, rho, Kernel);
  567. TERMR+=TR;
  568. TERMI+=TI;
  569. AA=BB;
  570. BB=BB+XINC;
  571. ++np;
  572. }
  573. }
  574. ++NumberPartialIntegrals;
  575. Sr(L) = TERMR;
  576. Si(L) = TERMI;
  577. Padecf(SUMR,SUMI,L);
  578. if ((std::abs(SUMR-LASTR) <= RERR*std::abs(SUMR)+AERR) &&
  579. (std::abs(SUMI-LASTI) <= RERR*std::abs(SUMI)+AERR)) {
  580. BESR=XSUMR+SUMR;
  581. BESI=XSUMI+SUMI;
  582. IERR=0;
  583. nps = std::max(np-1, nps);
  584. return;
  585. } else {
  586. LASTR=SUMR;
  587. LASTI=SUMI;
  588. ++NPB;
  589. A=B;
  590. B=ZeroJ(NPB,IORDER)/rho;
  591. ++L;
  592. }
  593. }
  594. BESR=XSUMR+SUMR;
  595. BESI=XSUMI+SUMI;
  596. IERR=1;
  597. nps = std::max(np, nps);
  598. return;
  599. }
  600. // GUSSIAN QUADRATURE ONLY IN THE INTERVAL IN WHICH K_0 IS LESS THAN
  601. // LAMBDA
  602. //std::cout << "NCNTRL " << NCNTRL << " NPCS " << NPCS << std::endl;
  603. for (int N=1; N<=NSTOP; ++N) {
  604. //std::cout << "NSTOP" << NSTOP << std::endl;
  605. if (NPCS == 1) {
  606. if ((NW==2) && (np > NPO)) {
  607. NW=1;
  608. }
  609. if (np > NTERM) {
  610. NW=0;
  611. }
  612. // Problem here, A not getting initialised
  613. //std::cout << "A " << A << " B " << B << std::endl;
  614. Besqud(A, B, TERMR, TERMI, NG, NW, IORDER, rho, Kernel);
  615. ++np;
  616. } else {
  617. TERMR=0.;
  618. TERMI=0.;
  619. Real XINC=(B-A)/NPCS;
  620. Real AA=A;
  621. Real BB=A+XINC;
  622. Real TR(0), TI(0);
  623. for (int I=1; I<NPCS; ++I) {
  624. if ((NW == 2) && (np > NPO)) {
  625. NW=1;
  626. }
  627. if (np > NTERM) {
  628. NW=0;
  629. }
  630. //std::cout << "AA " << AA << " BB " << BB << std::endl;
  631. Besqud(AA, BB, TR, TI, NG, NW, IORDER, rho, Kernel);
  632. TERMR+=TR;
  633. TERMI+=TI;
  634. AA=BB;
  635. BB=BB+XINC;
  636. ++np;
  637. }
  638. }
  639. if (B >= AORB) {
  640. BESR=XSUMR+TERMR;
  641. BESI=XSUMI+TERMI;
  642. IERR=0;
  643. nps = std::max(np, nps); // TI np is zero based
  644. return;
  645. } else {
  646. XSUMR = XSUMR + TERMR;
  647. XSUMI = XSUMI + TERMI;
  648. ++NPB;
  649. A = B;
  650. B = ZeroJ(NPB,IORDER)/rho;
  651. if (B >= AORB) B = AORB;
  652. }
  653. }
  654. return;
  655. }
  656. /////////////////////////////////////////////////////////////
  657. void GQChave::
  658. Besqud(const Real &LowerLimit, const Real &UpperLimit,
  659. Real &Besr, Real &Besi, const int &npoints,
  660. const int &NEW,
  661. const int &besselOrder, const Real &rho,
  662. KernelEM1DBase* Kernel) {
  663. const int NTERM = 100;
  664. // Filter Weights
  665. Eigen::Matrix<int , 7, 1> NWA;
  666. Eigen::Matrix<int , 7, 1> NWT;
  667. Eigen::Matrix<Real, 254, 1> FUNCT;
  668. Eigen::Matrix<Real, 64, 1> FR1;
  669. Eigen::Matrix<Real, 64, 1> FI1;
  670. Eigen::Matrix<Real, 64, 1> FR2;
  671. Eigen::Matrix<Real, 64, 1> FI2;
  672. Eigen::Matrix<Real, 64, 1> BES1;
  673. Eigen::Matrix<Real, 64, 1> BES2;
  674. // Init Vals
  675. NWT << 1, 3, 7, 15, 31, 63, 127;
  676. NWA << 1, 2, 4, 8, 16, 32, 64;
  677. // CHECK FOR TRIVIAL CASE
  678. if (LowerLimit >= UpperLimit) {
  679. Besr = 0.;
  680. Besi = 0.;
  681. //std::cout << "Lower limit > Upper limit " << NEW << std::endl;
  682. //std::cout << "Lower limit " << LowerLimit << "Upper Limit " << UpperLimit <<std::endl;
  683. return;
  684. }
  685. // Temp vars for function generation
  686. Real diff(0);
  687. Real FZEROR(0);
  688. Real FZEROI(0);
  689. Real SUM(0);
  690. int KN(0);
  691. int LA(0);
  692. int LK(0);
  693. int N(0);
  694. int K;
  695. int NW;
  696. Real ACUMR;
  697. Real ACUMI;
  698. Real BESF;
  699. switch (NEW) {
  700. // CONSTRUCT FUNCT FROM SAVED KERNELS
  701. case (2):
  702. diff = (UpperLimit-LowerLimit)/2.;
  703. FZEROR = kern(0, np);
  704. FZEROI = kern(1, np);
  705. K = std::min(Nk(np)-1, npoints-1);
  706. NW=NWT(K);
  707. LK=2;
  708. for (N=0; N<2*NW; N+=2) {
  709. FUNCT(N) = kern(LK , np) +
  710. kern(LK+2, np) ;
  711. FUNCT(N+1) = kern(LK+1, np) +
  712. kern(LK+3, np) ;
  713. LK=LK+4;
  714. }
  715. if (Nk(np) >= npoints) {
  716. ACUMR = _dot(NW, WT.tail(255-NW), 1, FUNCT, 2);
  717. ACUMI = _dot(NW, WT.tail(255-NW), 1, FUNCT.tail<253>(), 2);
  718. Besr = (ACUMR+WT(2*NW)*FZEROR)*diff;
  719. Besi = (ACUMI+WT(2*NW)*FZEROI)*diff;
  720. return;
  721. } else {
  722. // COMPUTE ADDITIONAL ORDERS BEFORE TAKING DOT PRODUCT
  723. SUM= (UpperLimit+LowerLimit)/2.;
  724. KN = K +1;
  725. LA = LK/2;
  726. }
  727. break;
  728. default:
  729. // SCALE FACTORS
  730. SUM = (UpperLimit+LowerLimit) / 2.;
  731. diff = (UpperLimit-LowerLimit) / 2.;
  732. // ONE POINT GAUSS
  733. //Kernel(SUM,FZEROR,FZEROI,TEM,GEOM,PARA);
  734. Complex tba = Kernel->BesselArg(SUM);
  735. //cout.precision(16);
  736. //std::cout << "SUM " << "\t" << SUM << "\t" << tba << "\n";
  737. ++NumberFunctionEvals;
  738. BESF = Jbess(SUM*rho,besselOrder);
  739. FZEROR = BESF*std::real(tba);//FZEROR;
  740. FZEROI = BESF*std::imag(tba);//FZEROI;
  741. if (NEW == 1) {
  742. karg(0, np) = SUM;
  743. kern(0, np) = FZEROR;
  744. kern(1, np) = FZEROI;
  745. LA=1;
  746. LK=2;
  747. }
  748. N = 0; // TI was 1, 0 based index
  749. KN = 0;
  750. } // End of new switch
  751. // STEP THROUGH GAUSS ORDERS (NG = Number Gauss)
  752. for (int K=KN; K<npoints; ++K) {
  753. // COMPUTE NEW FUNCTION VALUES
  754. int NA = NWA(K) - 1; // NWA assumes 1 based indexing
  755. //std::cout << "Besqud NA=" << NA << std::endl;
  756. for (int J=0; J<NWA(K); ++J) {
  757. Real X = WA(NA)*diff;
  758. Real SUMP = SUM+X;
  759. Real SUMM = SUM-X;
  760. ++NA;
  761. //std::cout << "Calling from here\n";
  762. //Kernel(SUMP, FR1(J), FI1(J), TEM, GEOM, PARA);
  763. //Kernel(SUMM, FR2(J), FI2(J), TEM, GEOM, PARA);
  764. Complex bes1 = Kernel->BesselArg(SUMP);
  765. //std::cout << "SUMP" << "\t" << SUMP << "\t" << bes1 << "\n";
  766. Complex bes2 = Kernel->BesselArg(SUMM);
  767. //std::cout << "SUMM" << "\t" << SUMM << "\t" << bes1 << "\n";
  768. FR1(J) = std::real(bes1);
  769. FI1(J) = std::imag(bes1);
  770. FR2(J) = std::real(bes2);
  771. FI2(J) = std::imag(bes2);
  772. NumberFunctionEvals +=2;
  773. BES1(J)=Jbess(SUMP*rho, besselOrder);
  774. BES2(J)=Jbess(SUMM*rho, besselOrder);
  775. if (NEW >= 1) {
  776. karg(LA , np) = SUMP;
  777. karg(LA+1, np) = SUMM;
  778. LA += 2;
  779. }
  780. }
  781. // COMPUTE PRODUCTS OF KERNELS AND BESSEL FUNCTIONS
  782. // TODO vectorize
  783. for (int J=0; J<NWA(K); ++J) {
  784. FR1(J ) = BES1(J)*FR1(J);
  785. FI1(J ) = BES1(J)*FI1(J);
  786. FR2(J ) = BES2(J)*FR2(J);
  787. FI2(J ) = BES2(J)*FI2(J);
  788. FUNCT(N ) = FR1(J)+FR2(J);
  789. FUNCT(N+1) = FI1(J)+FI2(J);
  790. N+=2;
  791. }
  792. if (NEW >= 1) {
  793. for (int J=0; J<NWA(K); ++J) {
  794. kern(LK , np) = FR1(J);
  795. kern(LK+1, np) = FI1(J);
  796. kern(LK+2, np) = FR2(J);
  797. kern(LK+3, np) = FI2(J);
  798. LK=LK+4;
  799. }
  800. }
  801. }
  802. // COMPUTE DOT PRODUCT OF WEIGHTS WITH INTEGRAND VALUES
  803. NW=NWT(npoints-1);
  804. ACUMR = _dot(NW, WT.tail(255-NW), 1, FUNCT , 2);
  805. ACUMI = _dot(NW, WT.tail(255-NW), 1, FUNCT.tail<253>(), 2);
  806. Besr = (ACUMR+WT(2*NW-1)*FZEROR)*diff;
  807. Besi = (ACUMI+WT(2*NW-1)*FZEROI)*diff;
  808. if (np <= NTERM) Nk(np) = npoints;
  809. return;
  810. }
  811. /////////////////////////////////////////////////////////////
  812. void GQChave::
  813. Padecf(Real &SUMR, Real &SUMI, const int &N) {
  814. if (N < 2) {
  815. // INITIALIZE FOR RECURSIVE CALCULATIONS
  816. //if (N == 0) {
  817. // Xr.setZero();
  818. // Xi.setZero();
  819. //Xr = Eigen::Matrix<Real, Eigen::Dynamic, 1>::Zero(100);
  820. //Xi = Eigen::Matrix<Real, Eigen::Dynamic, 1>::Zero(100);
  821. //}
  822. Dr(N) = Sr(N);
  823. Di(N) = Si(N);
  824. Real DENOM;
  825. switch (N) {
  826. case 0:
  827. DENOM = 1;
  828. Cfcor(N) = -(-1.*Dr(N)) / DENOM;
  829. Cfcoi(N) = -(-1.*Di(N)) / DENOM;
  830. break;
  831. default:
  832. DENOM = Dr(N-1)*Dr(N-1) +
  833. Di(N-1)*Di(N-1);
  834. Cfcor(N) = -(Dr(N-1)*Dr(N) +
  835. Di(N-1)*Di(N) ) / DENOM;
  836. Cfcoi(N) = -(Dr(N-1)*Di(N ) -
  837. Dr(N )*Di(N-1) ) / DENOM;
  838. }
  839. CF(SUMR, SUMI, Cfcor, Cfcoi, N);
  840. return;
  841. } else {
  842. int L = 2*(int)((N)/2) - 1;//+1;// - 1; // - 1; // - 1;
  843. // UPDATE X VECTORS FOR RECURSIVE CALCULATION OF CF COEFFICIENTS
  844. for (int K = L; K >= 3; K-= 2) {
  845. Xr(K) = Xr(K-1)+Cfcor(N-1)*Xr(K-2) -
  846. Cfcoi(N-1)*Xi(K-2);
  847. Xi(K) = Xi(K-1)+Cfcor(N-1)*Xi(K-2) +
  848. Cfcoi(N-1)*Xr(K-2);
  849. }
  850. Xr(1) = Xr(0)+Cfcor(N-1);
  851. Xi(1) = Xi(0)+Cfcoi(N-1);
  852. // INTERCHANGE ODD AND EVEN PARTS
  853. for (int K=0; K<L; K+=2) {
  854. Real T1 = Xr(K);
  855. Real T2 = Xi(K);
  856. Xr(K) = Xr(K+1);
  857. Xi(K) = Xi(K+1);
  858. Xr(K+1) = T1;
  859. Xi(K+1) = T2;
  860. }
  861. // COMPUTE FIRST COEFFICIENTS
  862. Dr(N) = Sr(N);
  863. Di(N) = Si(N);
  864. // Dr getting fucked up here
  865. for (int K=0; K<std::max(1, L/2+1); ++K) {
  866. Dr(N) += Sr(N-K-1)*Xr(2*K) -
  867. Si(N-K-1)*Xi(2*K);
  868. Di(N) += Si(N-K-1)*Xr(2*K) +
  869. Sr(N-K-1)*Xi(2*K);
  870. }
  871. // COMPUTE NEW CF COEFFICIENT
  872. Real DENOM = Dr(N-1)*Dr(N-1) +
  873. Di(N-1)*Di(N-1);
  874. //std::cout << "DENOM " << DENOM << std::endl;
  875. Cfcor(N)=-(Dr(N )*Dr(N-1) +
  876. Di(N)*Di(N-1))/DENOM;
  877. Cfcoi(N)=-(Dr(N-1)*Di(N ) -
  878. Dr(N)*Di(N-1))/DENOM;
  879. // EVALUATE CONTINUED FRACTION
  880. CF(SUMR,SUMI,Cfcor,Cfcoi,N);
  881. return;
  882. }
  883. }
  884. /////////////////////////////////////////////////////////////
  885. void GQChave::CF(Real& RESR, Real &RESI,
  886. Eigen::Matrix<Real, 100, 1> &CFCOR,
  887. Eigen::Matrix<Real, 100, 1> &CFCOI,
  888. const int &N) {
  889. ////////////////////////////////////////////////
  890. // ONE Seems sort of stupid, maybe use ? instead
  891. // TODO benchmark difference
  892. //RESR = ONE(N)+CFCOR(N);
  893. RESR = (!N ? 0:1) + CFCOR(N);
  894. RESI = CFCOI(N);
  895. for (int K=N-1; K>=0; --K) {
  896. Real DENOM=RESR*RESR + RESI*RESI;
  897. Real RESRO=RESR;
  898. //RESR=ONE(K)+(RESR*CFCOR(K)+RESI*CFCOI(K))/DENOM;
  899. RESR = (!K ? 0:1) + (RESR*CFCOR(K)+RESI*CFCOI(K))/DENOM;
  900. RESI=(RESRO*CFCOI(K)-RESI*CFCOR(K))/DENOM;
  901. }
  902. return;
  903. }
  904. /////////////////////////////////////////////////////////////
  905. Real GQChave::
  906. ZeroJ(const int &nzero, const int &besselOrder) {
  907. Real ZT1 = -1.e0/8.e0;
  908. Real ZT2 = 124.e0/1536.e0;
  909. Real ZT3 = -120928.e0/491520.e0;
  910. Real ZT4 = 401743168.e0/220200960.e0;
  911. Real OT1 = 3.e0/8.e0;
  912. Real OT2 = -36.e0/1536.e0;
  913. Real OT3 = 113184.e0/491520.e0;
  914. Real OT4 = -1951209.e0/220200960.e0;
  915. Real BETA(0);
  916. switch (besselOrder) {
  917. case 0:
  918. BETA=(nzero-.25e0)*PI;
  919. return BETA-ZT1/BETA-ZT2/std::pow(BETA,3) -
  920. ZT3/std::pow(BETA,5)-ZT4/std::pow(BETA,7);
  921. case 1:
  922. BETA=(nzero+.25e0)*PI;
  923. return BETA-OT1/BETA-OT2/std::pow(BETA,3) -
  924. OT3/std::pow(BETA,5)-OT4/std::pow(BETA,7);
  925. default:
  926. throw 77;
  927. }
  928. return 0.;
  929. }
  930. /////////////////////////////////////////////////////////////
  931. // Dot product allowing non 1 based incrementing
  932. Real GQChave::_dot(const int&n,
  933. const Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> &X1,
  934. const int &inc1,
  935. const Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> &X2,
  936. const int &inc2) {
  937. int k;
  938. if (inc2 > 0) {
  939. k=0;
  940. } else {
  941. k = n*std::abs(inc2);
  942. }
  943. Real dot=0.0;
  944. if (inc1 > 0) {
  945. for (int i=0; i<n; i += inc1) {
  946. dot += X1(i)*X2(k);
  947. k += inc2;
  948. }
  949. } else {
  950. for (int i=n-1; i>=0; i += inc1) {
  951. dot += X1(i)*X2(k);
  952. k += inc2;
  953. }
  954. }
  955. return dot;
  956. }
  957. /////////////////////////////////////////////////////////////
  958. //
  959. #ifdef HAVE_BOOST_SPECIAL_FUNCTIONS
  960. Real GQChave::Jbess(const Real &x, const int &IORDER) {
  961. switch (IORDER) {
  962. case 0:
  963. //return boost::math::detail::bessel_j0(X);
  964. return boost::math::cyl_bessel_j(0, x);
  965. break;
  966. case 1:
  967. //return boost::math::detail::bessel_j1(x);
  968. return boost::math::cyl_bessel_j(1, x);
  969. break;
  970. default:
  971. throw std::runtime_error("Non 0 or 1 Bessel argument specified in GQChave");
  972. }
  973. }
  974. #else
  975. Real GQChave::Jbess(const Real &x, const int &IORDER) {
  976. std::cerr << "GQChave requires boost special functions module";
  977. return 0;
  978. }
  979. #endif
  980. //////////////////////////////////////////////////////
  981. // Exception classes
  982. LowerGaussLimitGreaterThanUpperGaussLimit::
  983. LowerGaussLimitGreaterThanUpperGaussLimit() :
  984. runtime_error("LOWER GAUSS LIMIT > UPPER GAUSS LIMIT") {}
  985. } // ----- end of Lemma name -----