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.

FHTKey101.cpp 16KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274
  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/11/2014 03:42:53 PM
  11. * @version $Id$
  12. * @author Trevor Irons (ti)
  13. * @email Trevor.Irons@lemmasoftware.org
  14. * @copyright Copyright (c) 2014, Trevor Irons
  15. */
  16. #include "FHTKey101.h"
  17. namespace Lemma {
  18. // ==================== FRIEND METHODS =====================
  19. std::ostream &operator<<(std::ostream &stream, const FHTKey101 &ob) {
  20. stream << ob.Serialize() << "\n";
  21. return stream;
  22. }
  23. // ==================== STATIC CONST MEMBERS ============
  24. const Eigen::Matrix<Real, 101, 3> FHTKey101::WT101 =
  25. ( Eigen::Matrix<Real, 101, 3>() <<
  26. // Base J0 J1
  27. 5.5308437014783363e-04, 5.1818808036862153e-02, 4.1746363961646286e-06,
  28. 6.4259236035555788e-04, -2.9258492652779172e-01, -2.5993370843637242e-05,
  29. 7.4658580837667996e-04, 8.7289913080159265e-01, 9.2506838353721647e-05,
  30. 8.6740895730700251e-04, -1.8525183961240297e+00, -2.4318740269846489e-04,
  31. 1.0077854290485113e-03, 3.2096055660873608e+00, 5.2276050788732404e-04,
  32. 1.1708796207911744e-03, -4.9248158048315069e+00, -9.5803859755432658e-04,
  33. 1.3603680375478939e-03, 7.0219167445719828e+00, 1.5308357107396591e-03,
  34. 1.5805221687362171e-03, -9.4781099645811544e+00, -2.1584605764249656e-03,
  35. 1.8363047770289071e-03, 1.2098564895534539e+01, 2.7161344302423967e-03,
  36. 2.1334817700377098e-03, -1.4476661837664896e+01, -3.0845361397036195e-03,
  37. 2.4787521766663585e-03, 1.6146079632847474e+01, 3.2133980102756368e-03,
  38. 2.8798991580882430e-03, -1.6810370923047440e+01, -3.1305779240645847e-03,
  39. 3.3459654574712720e-03, 1.6506202923977266e+01, 2.9239120736637542e-03,
  40. 3.8874572434761303e-03, -1.5534747212675301e+01, -2.6680602560576316e-03,
  41. 4.5165809426126703e-03, 1.4267433361428399e+01, 2.4131379039476304e-03,
  42. 5.2475183991813838e-03, -1.2939053370802849e+01, -2.1549566068119440e-03,
  43. 6.0967465655156379e-03, 1.1628520257389562e+01, 1.9058200916563901e-03,
  44. 7.0834089290521185e-03, -1.0325503688577493e+01, -1.6533122444456158e-03,
  45. 8.2297470490200302e-03, 9.0574401975545840e+00, 1.4391633836110848e-03,
  46. 9.5616019305435132e-03, -7.8759096646654951e+00, -1.2363666069459527e-03,
  47. 1.1108996538242306e-02, 6.8446506699312666e+00, 1.0907227256323094e-03,
  48. 1.2906812580479873e-02, -5.9542991234257840e+00, -9.2712769282234358e-04,
  49. 1.4995576820477703e-02, 5.1830565431894602e+00, 8.2888526218099687e-04,
  50. 1.7422374639493515e-02, -4.4805481279340267e+00, -6.7675613963809832e-04,
  51. 2.0241911445804391e-02, 3.8608298760715378e+00, 6.3558641899842260e-04,
  52. 2.3517745856009110e-02, -3.3125736597743227e+00, -4.8591507683128695e-04,
  53. 2.7323722447292573e-02, 2.8648088091590589e+00, 5.1892328126100203e-04,
  54. 3.1745636378067953e-02, -2.4710577422247479e+00, -3.3190335155472070e-04,
  55. 3.6883167401240015e-02, 2.1450315408533416e+00, 4.6189809360776176e-04,
  56. 4.2852126867040180e-02, -1.8334353019592327e+00, -1.8326476919948130e-04,
  57. 4.9787068367863944e-02, 1.5850020809005190e+00, 4.7669939315323097e-04,
  58. 5.7844320874838456e-02, -1.3468969891769780e+00, -1.7380345519175726e-05,
  59. 6.7205512739749784e-02, 1.1802299801807625e+00, 5.9009214750893949e-04,
  60. 7.8081666001153169e-02, -9.9988174688478793e-01, 2.1645139320505843e-04,
  61. 9.0717953289412512e-02, 8.8503813256539254e-01, 8.5203135609860704e-04,
  62. 1.0539922456186433e-01, -7.3165077099639619e-01, 5.9560990889710838e-04,
  63. 1.2245642825298191e-01, 6.5909386165012596e-01, 1.3651107378683781e-03,
  64. 1.4227407158651359e-01, -5.3038429622108840e-01, 1.2606987195022564e-03,
  65. 1.6529888822158656e-01, 5.0119791866928076e-01, 2.3131990435032663e-03,
  66. 1.9204990862075413e-01, -3.8269840749371159e-01, 2.4572002938857678e-03,
  67. 2.2313016014842982e-01, 3.8622329595139565e-01, 4.0408430887840449e-03,
  68. 2.5924026064589156e-01, -2.6413330193320073e-01, 4.6278078809669602e-03,
  69. 3.0119421191220208e-01, 3.0353632730294311e-01, 7.1437051276295131e-03,
  70. 3.4993774911115533e-01, -1.7255655176111537e-01, 8.5723466849983495e-03,
  71. 4.0656965974059917e-01, 2.5119774955173280e-01, 1.2684763573137947e-02,
  72. 4.7236655274101469e-01, -9.8536724013140092e-02, 1.5643387012124067e-02,
  73. 5.4881163609402650e-01, 2.1755876398162183e-01, 2.2480928663078600e-02,
  74. 6.3762815162177333e-01, -3.4268716744688760e-02, 2.8124131719595692e-02,
  75. 7.4081822068171788e-01, 1.9925218957836366e-01, 3.9373999883133115e-02,
  76. 8.6070797642505781e-01, 1.8319735495994052e-02, 4.9445586323210833e-02,
  77. 1.0000000000000000e+00, 1.8963099219918259e-01, 6.7326587032714463e-02,
  78. 1.1618342427282831e+00, 5.7228135568732394e-02, 8.3541701874077570e-02,
  79. 1.3498588075760032e+00, 1.7266316834158107e-01, 1.0961776964853986e-01,
  80. 1.5683121854901687e+00, 6.8141322181604499e-02, 1.3118216481048459e-01,
  81. 1.8221188003905089e+00, 1.2468074033081850e-01, 1.6144181878027153e-01,
  82. 2.1170000166126748e+00, 2.1678852838804014e-02, 1.7704476632134886e-01,
  83. 2.4596031111569494e+00, 1.1114005126278175e-02, 1.9032654055463888e-01,
  84. 2.8576511180631639e+00, -1.0424806190544497e-01, 1.6333927603825021e-01,
  85. 3.3201169227365472e+00, -1.6473681128667569e-01, 1.1034442879456068e-01,
  86. 3.8574255306969740e+00, -2.3603153458714107e-01, -1.1282682226217962e-02,
  87. 4.4816890703380645e+00, -2.2093965064386309e-01, -1.4579378661401318e-01,
  88. 5.2069798271798486e+00, -7.5633542915308175e-02, -2.7560051939484942e-01,
  89. 6.0496474644129448e+00, 1.3358180223826097e-01, -2.3333751023296440e-01,
  90. 7.0286875805892928e+00, 3.3949946683489907e-01, -5.9615109981928909e-03,
  91. 8.1661699125676499e+00, 1.2815187347111337e-01, 3.2380356528743487e-01,
  92. 9.4877358363585262e+00, -2.2988301163130145e-01, 2.1970582679149839e-01,
  93. 1.1023176380641601e+01, -3.3425750283830036e-01, -2.8302631898505476e-01,
  94. 1.2807103782663029e+01, 4.4350537603751145e-01, -2.2706962776315812e-01,
  95. 1.4879731724872830e+01, -8.8943820418669634e-02, 4.7351258711649669e-01,
  96. 1.7287781840567639e+01, -1.8646079800862164e-01, -3.5417866045339608e-01,
  97. 2.0085536923187668e+01, 2.4494657419883456e-01, 1.6582378621661517e-01,
  98. 2.3336064580942711e+01, -2.0169045144304532e-01, -5.1221107687420810e-02,
  99. 2.7112638920657883e+01, 1.4661206655803793e-01, 3.5484710794282063e-03,
  100. 3.1500392308747923e+01, -1.0564640256629308e-01, 1.1614894738175324e-02,
  101. 3.6598234443677974e+01, 7.8487036158411266e-02, -1.4776256489774963e-02,
  102. 4.2521082000062783e+01, -6.0533748554794195e-02, 1.4194293010344678e-02,
  103. 4.9402449105530167e+01, 4.8243702864613026e-02, -1.2609914713728791e-02,
  104. 5.7397457045446188e+01, -3.9347085461672239e-02, 1.0854082569077052e-02,
  105. 6.6686331040925154e+01, 3.2489421068688185e-02, -9.2082270347668624e-03,
  106. 7.7478462925260828e+01, -2.6982535264609676e-02, 7.7880648968821110e-03,
  107. 9.0017131300521811e+01, 2.2530196869950080e-02, -6.6089268276601255e-03,
  108. 1.0458498557711414e+02, -1.8947601435977172e-02, 5.6159700146012417e-03,
  109. 1.2151041751873485e+02, 1.6017457455233256e-02, -4.7432063636330263e-03,
  110. 1.4117496392147686e+02, -1.3519422964231352e-02, 3.9609893136500528e-03,
  111. 1.6402190729990167e+02, 1.1318328695851308e-02, -3.2753768787828126e-03,
  112. 1.9056626845862999e+02, -9.3839154851987506e-03, 2.6945716306237466e-03,
  113. 2.2140641620418697e+02, 7.7256174599644543e-03, -2.2042319512562145e-03,
  114. 2.5723755590577474e+02, -6.3199190181430056e-03, 1.7733135127602079e-03,
  115. 2.9886740096706029e+02, 5.0986473606509210e-03, -1.3773733375919272e-03,
  116. 3.4723438047873447e+02, -3.9929354358609528e-03, 1.0137885881408651e-03,
  117. 4.0342879349273511e+02, 2.9783618569305080e-03, -6.9751488719640639e-04,
  118. 4.6871738678241655e+02, -2.0802747552777045e-03, 4.4545125438205673e-04,
  119. 5.4457191012592898e+02, 1.3429398653973710e-03, -2.6373575275283686e-04,
  120. 6.3270229281225352e+02, -7.9345560431147465e-04, 1.4519595835116970e-04,
  121. 7.3509518924197266e+02, 4.2530878973518046e-04, -7.4640168948626184e-05,
  122. 8.5405876252615155e+02, -2.0457624610548748e-04, 3.5826939350044606e-05,
  123. 9.9227471560502534e+02, 8.6750842046621606e-05, -1.5854604622793124e-05,
  124. 1.1528587427833875e+03, -3.1406113278782205e-05, 6.2412989356183992e-06,
  125. 1.3394307643944169e+03, 9.1385289919391839e-06, -2.0349222128793154e-06,
  126. 1.5561965278371533e+03, -1.8949818224609619e-06, 4.8072849734177625e-07,
  127. 1.8080424144560632e+03, 2.0794387557779629e-07, -6.0462736574031818e-08).finished();
  128. // ==================== LIFECYCLE =======================
  129. //--------------------------------------------------------------------------------------
  130. // Class: FHTKey101
  131. // Method: FHTKey101
  132. // Description: constructor (locked)
  133. //--------------------------------------------------------------------------------------
  134. FHTKey101::FHTKey101 ( const ctor_key& key ) : HankelTransform( key ) {
  135. } // ----- end of method FHTKey101::FHTKey101 (constructor) -----
  136. //--------------------------------------------------------------------------------------
  137. // Class: FHTKey101
  138. // Method: FHTKey101
  139. // Description: constructor (protected)
  140. //--------------------------------------------------------------------------------------
  141. FHTKey101::FHTKey101( const YAML::Node& node, const ctor_key& key ) : HankelTransform(node, key) {
  142. }
  143. //--------------------------------------------------------------------------------------
  144. // Class: FHTKey101
  145. // Method: New()
  146. // Description: public constructor
  147. //--------------------------------------------------------------------------------------
  148. std::shared_ptr<FHTKey101> FHTKey101::NewSP() {
  149. return std::make_shared< FHTKey101 >( ctor_key() );
  150. }
  151. //--------------------------------------------------------------------------------------
  152. // Class: FHTKey101
  153. // Method: ~FHTKey101
  154. // Description: destructor (protected)
  155. //--------------------------------------------------------------------------------------
  156. FHTKey101::~FHTKey101 () {
  157. } // ----- end of method FHTKey101::~FHTKey101 (destructor) -----
  158. //--------------------------------------------------------------------------------------
  159. // Class: FHTKey101
  160. // Method: DeSerialize
  161. // Description: Factory method, converts YAML node into object
  162. //--------------------------------------------------------------------------------------
  163. std::shared_ptr<FHTKey101> FHTKey101::DeSerialize( const YAML::Node& node ) {
  164. if (node.Tag() != "FHTKey101") {
  165. throw DeSerializeTypeMismatch( "FHTKey101", node.Tag());
  166. }
  167. return std::make_shared<FHTKey101> ( node, ctor_key() );
  168. }
  169. //--------------------------------------------------------------------------------------
  170. // Class: FHTKey101
  171. // Method: Serialize
  172. // Description: Converts object into Serialized version
  173. //--------------------------------------------------------------------------------------
  174. YAML::Node FHTKey101::Serialize() const {
  175. YAML::Node node = HankelTransform::Serialize();
  176. node.SetTag( GetName() );
  177. //node["LayerConductivity"] = LayerConductivity;
  178. return node;
  179. }
  180. //--------------------------------------------------------------------------------------
  181. // Class: FHTKey101
  182. // Method: GetName
  183. // Description: Class identifier
  184. //--------------------------------------------------------------------------------------
  185. inline std::string FHTKey101::GetName ( ) const {
  186. return CName;
  187. } // ----- end of method FHTKey101::GetName -----
  188. //--------------------------------------------------------------------------------------
  189. // Class: FHTKey101
  190. // Method: Zgauss
  191. //--------------------------------------------------------------------------------------
  192. Complex FHTKey101::Zgauss ( const int &ikk, const EMMODE &imode,
  193. const int &itype, const Real &rho,
  194. const Real &wavef, KernelEM1DBase* Kernel ) {
  195. return Zans(0, Kernel->GetManagerIndex());
  196. } // ----- end of method FHTKey101::ComputeRelated -----
  197. //--------------------------------------------------------------------------------------
  198. // Class: FHTKey101
  199. // Method: ComputeRelated
  200. //--------------------------------------------------------------------------------------
  201. void FHTKey101::ComputeRelated ( const Real& rho, std::shared_ptr<KernelEM1DBase> Kernel ) {
  202. return ;
  203. } // ----- end of method FHTKey101::ComputeRelated -----
  204. //--------------------------------------------------------------------------------------
  205. // Class: FHTKey101
  206. // Method: ComputeRelated
  207. //--------------------------------------------------------------------------------------
  208. void FHTKey101::ComputeRelated ( const Real& rho, std::vector< std::shared_ptr<KernelEM1DBase> > KernelVec ) {
  209. return ;
  210. } // ----- end of method FHTKey101::ComputeRelated -----
  211. //--------------------------------------------------------------------------------------
  212. // Class: FHTKey101
  213. // Method: ComputeRelated
  214. //--------------------------------------------------------------------------------------
  215. void FHTKey101::ComputeRelated ( const Real& rho, std::shared_ptr<KernelEM1DManager> KernelManager ) {
  216. //kernelVec = KernelManager->GetSTLVector();
  217. int nrel = (int)(KernelManager->GetSTLVector().size());
  218. Eigen::Matrix<Complex, 101, Eigen::Dynamic > Zwork;
  219. // TODO, if we want to allow lagged, then 1 below should be nlag
  220. Zans= Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(1, nrel);
  221. Zwork.resize(101, nrel);
  222. VectorXr lambda = WT101.col(0).array()/rho;
  223. int NumFun = 0;
  224. int idx = 0;
  225. // Get Kernel values
  226. for (int ir=0; ir<lambda.size(); ++ir) {
  227. // irelated loop
  228. ++NumFun;
  229. KernelManager->ComputeReflectionCoeffs(lambda(ir), idx, rho);
  230. for (int ir2=0; ir2<nrel; ++ir2) {
  231. // Zwork* needed due to sign convention of filter weights
  232. Zwork(ir, ir2) = std::conj(KernelManager->GetSTLVector()[ir2]->RelBesselArg(lambda(ir)));
  233. }
  234. }
  235. // We diverge slightly from Key here, each kernel is evaluated seperately, whereby instead
  236. // they prefer to sum them. The reason is that all those terms have been removed from the kernels
  237. // in the interests of making them as generic and reusable as possible. This approach requires slightly
  238. // more multiplies, but the same number of kernel evaluations, which is the expensive part.
  239. // Inner product and scale
  240. for (int ir2=0; ir2<nrel; ++ir2) {
  241. Zans(0, ir2) = Zwork.col(ir2).dot(WT101.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder() + 1))/rho;
  242. }
  243. return ;
  244. } // ----- end of method FHTKey101::ComputeRelated -----
  245. } // ----- end of Lemma name -----