Lemma is an Electromagnetics API
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.

InstrumentTEM.cpp 16KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449
  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 M. Andy Kass, Edited Trevor Irons 02/28/2014
  8. @date 02/10/2011
  9. @version $Id: instrumenttem.cpp 266 2015-04-01 03:24:00Z tirons $
  10. **/
  11. #include "InstrumentTEM.h"
  12. namespace Lemma {
  13. std::ostream &operator<<(std::ostream &stream,
  14. const InstrumentTem &ob) {
  15. stream << *(LemmaObject*)(&ob);
  16. return stream;
  17. }
  18. // ======================= Lifecycle ==================================
  19. InstrumentTem::InstrumentTem(const std::string &name) : Instrument(name),
  20. Antenna(NULL), Receiver(NULL), EarthMod(NULL), Dipole(NULL), bipolarWaveform(true),
  21. NPreviousPulse(0), RefTime(0.0), PulseT(.06), TxFreq(0.30), RType(INDUCTIVE) {
  22. }
  23. InstrumentTem::~InstrumentTem() {
  24. if (NumberOfReferences != 0) {
  25. throw DeleteObjectWithReferences(this);
  26. }
  27. if (this->Antenna != NULL) {
  28. this->Antenna->DetachFrom(this);
  29. }
  30. if (this->Dipole != NULL) {
  31. this->Dipole->DetachFrom(this);
  32. }
  33. if (this->EarthMod != NULL) {
  34. this->EarthMod->DetachFrom(this);
  35. }
  36. if (this->Receiver != NULL) {
  37. this->Receiver->DetachFrom(this);
  38. }
  39. }
  40. InstrumentTem* InstrumentTem::New() {
  41. InstrumentTem* Obj = new InstrumentTem("InstrumentTem");
  42. Obj->AttachTo(Obj);
  43. return Obj;
  44. }
  45. void InstrumentTem::Delete() {
  46. this->DetachFrom(this);
  47. }
  48. void InstrumentTem::Release() {
  49. delete this;
  50. }
  51. // ======================= Operations =================================
  52. void InstrumentTem::MakeDirectCalculation( const HANKELTRANSFORMTYPE& hType) {
  53. //int temp;
  54. int NumTimes;
  55. Real answer;
  56. EMEarth1D *EMEarthLocal=EMEarth1D::New();
  57. EMEarthLocal->AttachLayeredEarthEM(this->EarthMod);
  58. //DigitalFilterCosTrans *DFIntegrate=DigitalFilterCosTrans::New();
  59. DigitalFilterSinTrans *DFIntegrate=DigitalFilterSinTrans::New();
  60. TemIntegrationKernel *TEMIntKern=TemIntegrationKernel::New();
  61. TEMIntKern->SetEMEarth1D(EMEarthLocal);
  62. //TEMIntKern->SetTransmitter(this->Antenna);
  63. TEMIntKern->SetReceiver(this->Receiver);
  64. DFIntegrate->AttachKernel(TEMIntKern);
  65. if (this->Dipole != NULL) { //AND this->Antenna == NULL
  66. EMEarthLocal->AttachDipoleSource(this->Dipole);
  67. TEMIntKern->SetDipole(this->Dipole);
  68. }
  69. if (this->Antenna != NULL) { //AND this->Dipole == NULL
  70. EMEarthLocal->AttachWireAntenna(this->Antenna);
  71. TEMIntKern->SetTransmitter(this->Antenna);
  72. }
  73. EMEarthLocal->SetFieldsToCalculate(H);
  74. //EMEarthLocal->SetHankelTransformMethod(DIGITALFILTERING);
  75. //EMEarthLocal->SetHankelTransformMethod(FHTKEY201);
  76. EMEarthLocal->SetHankelTransformMethod(hType);
  77. //EMEarthLocal->SetHankelTransformMethod(GAUSSIANQUADRATURE);
  78. //EMEarthLocal->AttachReceiverPoints(this->Receiver);
  79. // For testing with one time only
  80. //DFIntegrate->Compute(1,1,1e-16);
  81. //std::cout<<DFIntegrate->GetAnswer()<<std::endl;
  82. //
  83. NumTimes = this->TimeGates.size();
  84. //Testing
  85. //std::cout << "NumTimes: " << NumTimes << std::endl;
  86. this->ModelledData.resize(NumTimes,2);
  87. for (int ii=0;ii<NumTimes;ii++) {
  88. //TESTING!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  89. //RELAXING THE VALUE!
  90. DFIntegrate->Compute(this->TimeGates(ii), 1, 1e-7);
  91. //Testing
  92. //std::cout << "DFIntegrate done computing" << std::endl;
  93. answer = DFIntegrate->GetAnswer()(0);
  94. answer = (2./PI)*answer;
  95. this->ModelledData(ii,0) = this->TimeGates(ii);
  96. this->ModelledData(ii,1) = answer;
  97. //std::cout << " " << TimeGates(ii) << " " <<
  98. // answer << std::endl;
  99. //std::cout << " " << this->ModelledData(ii,0) << " " <<
  100. // this->ModelledData(ii,1) << std::endl;
  101. // Call cosine transform which calls argument
  102. }
  103. //
  104. TEMIntKern->Delete();
  105. EMEarthLocal->Delete();
  106. DFIntegrate->Delete();
  107. }
  108. void InstrumentTem::MakeLaggedCalculation(const HANKELTRANSFORMTYPE& hType) {
  109. EMEarth1D *EMEarthLocal=EMEarth1D::New();
  110. EMEarthLocal->AttachLayeredEarthEM(this->EarthMod);
  111. EMEarthLocal->SetFieldsToCalculate(H);
  112. EMEarthLocal->SetHankelTransformMethod(hType);
  113. EMEarthLocal->AttachReceiverPoints(this->Receiver);
  114. TemIntegrationKernel *TEMIntKern=TemIntegrationKernel::New();
  115. TEMIntKern->SetReceiver(this->Receiver);
  116. TEMIntKern->SetEMEarth1D(EMEarthLocal);
  117. if (this->Dipole != NULL) { //AND this->Antenna == NULL
  118. EMEarthLocal->AttachDipoleSource(this->Dipole);
  119. TEMIntKern->SetDipole(this->Dipole);
  120. Dipole->SetNumberOfFrequencies(1);
  121. }
  122. if (this->Antenna != NULL) { //AND this->Dipole == NULL
  123. EMEarthLocal->AttachWireAntenna(this->Antenna);
  124. Antenna->SetCurrent(1);
  125. TEMIntKern->SetTransmitter(this->Antenna);
  126. Antenna->SetNumberOfFrequencies(1);
  127. }
  128. // TODO IMPORTANT!!
  129. // Compute number of lagged convolutions to do, TODO, need to make sure we do these
  130. // TO THE LATEST time needed for multiple pulse subtraction!!!
  131. int ilag = 1;
  132. Real sub = 2.1*PulseT + 2.*NPreviousPulse*PulseT;
  133. if (!bipolarWaveform) sub = NPreviousPulse*PulseT;
  134. //TimeGates.array() += RefTime;
  135. //RefTime = 0;
  136. // This is a critical line, for reasons I don't understand, it's important to calculate past the
  137. // theoretical GateWidths(0)/2.
  138. while ( (TimeGates.tail(1)(0)+sub)*std::pow(1./exp(.1), ilag-1) > TimeGates(0)-GateWidths(0) ) ++ ilag;
  139. //DigitalFilterSinTrans *DFIntegrate=DigitalFilterSinTrans::New(); // use with step response kernel
  140. DigitalFilterCosTrans *DFIntegrate=DigitalFilterCosTrans::New(); // use with impulse response kernel
  141. DFIntegrate->AttachKernel(TEMIntKern);
  142. DFIntegrate->SetNumConv(ilag);
  143. DFIntegrate->Compute( this->TimeGates.tail(1)(0)+sub, 1, 1e-7 ); // rho ntol, tol=1e-7 to 1e-13
  144. CubicSplineInterpolator *Spline = CubicSplineInterpolator::New();
  145. //Spline->SetKnots( DFIntegrate->GetAbscissaArguments().array(), (2./PI) * DFIntegrate->GetAnswer().col(0).array() );
  146. Spline->SetKnots( DFIntegrate->GetAbscissaArguments().array(), (.25/PI) * DFIntegrate->GetAnswer().col(0).array() );
  147. this->ModelledData.resize(TimeGates.size(), 2);
  148. this->ModelledData.col(0) = TimeGates;
  149. this->ModelledData.col(1).setZero();
  150. // This is non-extrapolated times and values
  151. //this->ModelledData.resize(ilag, 2);
  152. // this->ModelledData.col(0) = DFIntegrate->GetAbscissaArguments(); // this->TimeGates(ii);
  153. // this->ModelledData.col(1) = (2./PI)*DFIntegrate->GetAnswer();
  154. //std::cout << (2./PI)*DFIntegrate->GetAnswer().transpose() << std::endl;
  155. // Take into account receiver gate integration and tx waveform
  156. FoldAndConvolve(Spline);
  157. // Direct caclulation, step response to use need re reset temintegrationkernel back to step!
  158. // this->ModelledData.col(1) = (2./PI)*Spline->InterpolateOrderedSet( TimeGates );
  159. Spline->Delete();
  160. TEMIntKern->Delete();
  161. EMEarthLocal->Delete();
  162. DFIntegrate->Delete();
  163. }
  164. // ======================= Access =====================================
  165. void InstrumentTem::SetPulse( const VectorXr& Amp, const VectorXr& times, bool bipolar ) {
  166. bipolarWaveform = bipolar;
  167. assert(Amp.size() == times.size());
  168. TxAmp = Amp;
  169. TxAbs = times;
  170. TxDiff = VectorXr::Zero(TxAmp.size());
  171. TxDelta = VectorXr::Zero(TxAmp.size());
  172. for (int ip=1; ip<TxAmp.size(); ++ ip) {
  173. TxDiff(ip) = (TxAmp(ip-1)-TxAmp(ip )) /
  174. (TxAbs(ip )-TxAbs(ip-1)) ;
  175. TxDelta(ip) = (TxAmp(ip-1)-TxAmp(ip )) ;
  176. }
  177. /*
  178. std::cout << "diff " << TxDiff.transpose() << std::endl;
  179. std::cout << "delta " << TxDelta.transpose() << std::endl;
  180. std::cout << "amp " << TxAmp.transpose() << std::endl;
  181. std::cout << "abs " << TxAbs.transpose() << std::endl;
  182. */
  183. }
  184. void InstrumentTem::EMEarthModel(LayeredEarthEM* Earth) {
  185. if (this->EarthMod != NULL) {
  186. this->EarthMod->DetachFrom(this);
  187. }
  188. //std::cout << "InstrumentTem: EarthMod null test" << std::endl;
  189. Earth->AttachTo(this);
  190. //std::cout << "InstrumentTem: EarthMod attached" << std::endl;
  191. this->EarthMod = Earth;
  192. //std::cout << "InstrumentTem: EMEarthModel Attached" << std::endl;
  193. }
  194. void InstrumentTem::SetTransmitLoop(WireAntenna *antennae) {
  195. if (this->Antenna != NULL) {
  196. this->Antenna->DetachFrom(this);
  197. }
  198. antennae->AttachTo(this);
  199. this->Antenna = antennae;
  200. if (this->Dipole != NULL) {
  201. this->Dipole->DetachFrom(this);
  202. }
  203. //std::cout << "InstrumentTem: TransmitLoop Attached" << std::endl;
  204. }
  205. void InstrumentTem::SetDipoleSource(DipoleSource* dipolesource) {
  206. if (this->Dipole != NULL) {
  207. this->Dipole->DetachFrom(this);
  208. }
  209. dipolesource->AttachTo(this);
  210. this->Dipole = dipolesource;
  211. if (this->Antenna != NULL) {
  212. this->Antenna->DetachFrom(this);
  213. }
  214. }
  215. void InstrumentTem::SetReceiver(ReceiverPoints* Receivers) {
  216. if (this->Receiver != NULL) {
  217. this->Receiver->DetachFrom(this);
  218. }
  219. Receivers->AttachTo(this);
  220. this->Receiver = Receivers;
  221. //std::cout << "InstrumentTem: Receivers Attached" << std::endl;
  222. }
  223. void InstrumentTem::SetTimes(const VectorXr &times) {
  224. this->TimeGates = times;
  225. }
  226. //--------------------------------------------------------------------------------------
  227. // Class: InstrumentTem
  228. // Method: SetReceiverType
  229. //--------------------------------------------------------------------------------------
  230. void InstrumentTem::SetReceiverType ( const ReceiverType& rt ) {
  231. RType = rt;
  232. return ;
  233. } // ----- end of method InstrumentTem::SetReceiverType -----
  234. //--------------------------------------------------------------------------------------
  235. // Class: InstrumentTEM
  236. // Method: SetTimeGates
  237. //--------------------------------------------------------------------------------------
  238. void InstrumentTem::SetTimeGates ( const VectorXr &centres, const VectorXr& widths ) {
  239. assert (centres.size() == widths.size());
  240. TimeGates = centres;
  241. GateWidths = widths;
  242. return ;
  243. } // ----- end of method InstrumentTEM::SetTimeGates -----
  244. int InstrumentTem::GetNumberOfTimes() {
  245. return this->TimeGates.size();
  246. }
  247. MatrixXr InstrumentTem::GetMeasurements() {
  248. return this->ModelledData;
  249. }
  250. //--------------------------------------------------------------------------------------
  251. // Class: InstrumentTem
  252. // Method: SetReferenceTime
  253. //--------------------------------------------------------------------------------------
  254. void InstrumentTem::SetReferenceTime ( const Real& rtime ) {
  255. RefTime = rtime;
  256. return ;
  257. } // ----- end of method InstrumentTem::SetReferenceTime -----
  258. //--------------------------------------------------------------------------------------
  259. // Class: InstrumentTem
  260. // Method: FoldAndConvolve
  261. //--------------------------------------------------------------------------------------
  262. void InstrumentTem::FoldAndConvolve ( CubicSplineInterpolator* Spline ) {
  263. const static Real GLW[3] = {.5555556, .8888889, .5555556};
  264. const static Real GLX[3] = {-.7745967, 0., .7745967};
  265. // TODO BUG check against LEROI!!!
  266. //SubtractPrevious(Spline);
  267. switch (RType) {
  268. case INDUCTIVE:
  269. // Differentiate to compute impulse response for step input
  270. // Iterate over time gate
  271. for (int ig=0; ig<TimeGates.size(); ++ig) {
  272. this->ModelledData.col(1)[ig] =
  273. ( ( ConvolveWaveform(Spline, RefTime + TimeGates[ig] - GateWidths[ig]/2.) -
  274. ConvolveWaveform(Spline, RefTime + TimeGates[ig] + GateWidths[ig]/2.) ) / GateWidths[ig]);
  275. }
  276. break;
  277. case MAGNETOMETER:
  278. for (int ig=0; ig<TimeGates.size(); ++ig) {
  279. /*
  280. TC(2) = (TCLS(JT) + TOPN(JT)) / 2.
  281. TC(1) = TC(2) + HWIDTH * GLX(1)
  282. TC(3) = TC(2) + HWIDTH * GLX(3)
  283. */
  284. this->ModelledData.col(1)[ig] = 0;
  285. for (int ii=0; ii<3; ++ii) {
  286. //Real tt = RefTime + TimeGates[ig] + GLX[ii]*(GateWidths[ig]/2.);
  287. Real tt = RefTime + TimeGates[ig] + GLX[ii]*(GateWidths[ig]/2.);
  288. this->ModelledData.col(1)[ig] += ConvolveWaveform(Spline, tt) * GLW[ii]/2.;
  289. }
  290. }
  291. break;
  292. case ELECTRIC:
  293. throw std::runtime_error("InstrumentTEM->FoldAndConvolve with Electic Receivers is not implimented");
  294. break;
  295. }
  296. return ;
  297. } // ----- end of method InstrumentTem::FoldAndConvolve -----
  298. //--------------------------------------------------------------------------------------
  299. // Class: InstrumentTem
  300. // Method: ConvolveWaveform
  301. //--------------------------------------------------------------------------------------
  302. Real InstrumentTem::ConvolveWaveform ( CubicSplineInterpolator* Spline, const Real& t ) {
  303. static const Real T0MIN = 1e-7;
  304. Real cnv(0);
  305. Real tf = t - Spline->GetKnotAbscissa()[0];
  306. bool der = false;
  307. Real tend(0);
  308. for (int ip=1; ip<TxAmp.size(); ++ip) {
  309. if (TxAbs(ip) < T0MIN) continue;
  310. if (TxAbs(ip-1) > tf) break;
  311. Real tb = t - std::min(tf, TxAbs(ip));
  312. Real delt = TxAbs(ip) - TxAbs(ip-1);
  313. der = false;
  314. if (delt > T0MIN) {
  315. tend = t - TxAbs(ip-1);
  316. der = true;
  317. }
  318. if (der) {
  319. //std::cout.precision(16);
  320. // TODO Integrate seems to be kind of screwy!! Affecting early/on times
  321. cnv += TxDiff(ip) * Spline->Integrate(tb, tend, 41200); //CUBINT (TRP,YPLS,NTYPLS,TB,TEND)
  322. // TODO bug in direct spline integration, works sometimes
  323. //cnv += TxDiff(ip) * Spline->Integrate(tb, tend);
  324. //cnv += Spline->Integrate(tb, tend, 6264); //CUBINT (TRP,YPLS,NTYPLS,TB,TEND)
  325. //Real temp = Spline->Integrate(tb, tend, 10036); //CUBINT (TRP,YPLS,NTYPLS,TB,TEND)
  326. //Real temp2 = Spline->Integrate(tb, tend); //CUBINT (TRP,YPLS,NTYPLS,TB,TEND)
  327. //std::cout << "der " << TxDiff(ip) << "\t" << std::scientific << tb << "\t" << tend << "\t" << MU0*temp2 << "\t" << MU0*temp << "\n";
  328. //cnv += temp;
  329. } else {
  330. cnv += TxDelta(ip) * Spline->Interpolate(tb); //CUBVAL (TRP,YPLS,NTYPLS,TB)
  331. }
  332. }
  333. return cnv;
  334. } // ----- end of method InstrumentTem::ConvolveWaveform -----
  335. //--------------------------------------------------------------------------------------
  336. // Class: InstrumentTem
  337. // Method: SubtractPrevious
  338. //--------------------------------------------------------------------------------------
  339. void InstrumentTem::SubtractPrevious ( CubicSplineInterpolator *Spline ) {
  340. // Start by calculating impulse response of this pulse plus any arbitrary number of previous pulses at this
  341. // repetition frequency
  342. VectorXr aTimeGates = Spline->GetKnotAbscissa();
  343. Real pol (1.);
  344. if (bipolarWaveform) pol = -1.;
  345. VectorXr PrevPulse = aTimeGates.array() + PulseT;
  346. // TODO, if not bipolar what does convention say, is this added or should this be zero and handled only by nPreviousPulse?
  347. VectorXr Fold = Spline->GetKnotOrdinate( ); // +
  348. //VectorXr Fold = pol*Spline->InterpolateOrderedSet( PrevPulse ) ;
  349. VectorXr PrevPulse1; // = PrevPulse.array(); // + PulseT;
  350. for (int ip=0; ip<NPreviousPulse; ++ip) {
  351. PrevPulse1.array() += PrevPulse.array() + PulseT;
  352. PrevPulse.array() += PrevPulse1.array() + PulseT;
  353. Fold += Spline->InterpolateOrderedSet( PrevPulse1 ) + pol*Spline->InterpolateOrderedSet( PrevPulse ) ;
  354. }
  355. Spline->ResetKnotOrdinate(Fold);
  356. return ;
  357. } // ----- end of method InstrumentTem::SubtractPrevious -----
  358. } // namespace Lemma