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.

DipoleSource.cpp 17KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498
  1. // ===========================================================================
  2. //
  3. // Filename: utdipolesource.cpp
  4. //
  5. // Description:
  6. //
  7. // Version: 0.0
  8. // Created: 12/02/2009 11:57:14 AM
  9. // Revision: none
  10. // Compiler: g++ (c++)
  11. //
  12. // Author: Trevor Irons (ti)
  13. //
  14. // Organisation: Colorado School of Mines (CSM)
  15. // United States Geological Survey (USGS)
  16. //
  17. // Email: tirons@mines.edu, tirons@usgs.gov
  18. //
  19. // This program is free software: you can redistribute it and/or modify
  20. // it under the terms of the GNU General Public License as published by
  21. // the Free Software Foundation, either version 3 of the License, or
  22. // (at your option) any later version.
  23. //
  24. // This program is distributed in the hope that it will be useful,
  25. // but WITHOUT ANY WARRANTY; without even the implied warranty of
  26. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  27. // GNU General Public License for more details.
  28. //
  29. // You should have received a copy of the GNU General Public License
  30. // along with this program. If not, see <http://www.gnu.org/licenses/>.
  31. //
  32. // ===========================================================================
  33. #include <iostream>
  34. #include <fstream>
  35. #include "DipoleSource.h"
  36. // For testing purposes disable VTK and run scale.sh
  37. #undef LEMMAUSEVTK
  38. #undef KIHALEE_EM1D
  39. #ifdef LEMMAUSEVTK
  40. #include "vtkRenderer.h"
  41. #include "vtkRenderWindow.h"
  42. #include "vtkRenderWindowInteractor.h"
  43. #include "vtkDoubleArray.h"
  44. #include "vtkXYPlotActor.h"
  45. #include "vtkXYPlotWidget.h"
  46. #include "vtkProperty2D.h"
  47. #endif
  48. #include "timer.h"
  49. using namespace Lemma;
  50. int main() {
  51. //////////////////
  52. // TODO List of Dipole settings that aren't working
  53. // Do calculation
  54. jsw_timer timer;
  55. timer.begin();
  56. // Test with a single dipole
  57. DipoleSource *dipole = DipoleSource::New();
  58. //dipole->SetMoment(2);
  59. ///////////////////////////////////
  60. // E field Tests
  61. // Sems to be working!
  62. //dipole->SetType(GROUNDEDELECTRICDIPOLE);
  63. //dipole->SetPolarisation(XPOLARISATION);
  64. //dipole->SetPolarisation(YPOLARISATION);
  65. //dipole->SetPolarisation(ZPOLARISATION);
  66. //dipole->SetAzimuth(0); // X - Northing
  67. //dipole->SetAzimuth(90); // Y - Easting
  68. dipole->SetType(GROUNDEDELECTRICDIPOLE);
  69. //dipole->SetPolarisation(1., 0.0, 1.0);
  70. // dipole->SetPolarisation(XPOLARISATION);
  71. dipole->SetPolarisation(YPOLARISATION);
  72. //dipole->SetPolarisation(ZPOLARISATION);
  73. /////////////
  74. //dipole->SetType(MAGNETICDIPOLE);
  75. //dipole->SetPolarisation(0., 0.0, 1.0);
  76. //dipole->SetPolarisation(XPOLARISATION);
  77. //dipole->SetPolarisation(YPOLARISATION);
  78. //dipole->SetPolarisation(ZPOLARISATION);
  79. //dipole->SetMoment(1);
  80. //dipole->SetLocation(1,1,-.0100328);
  81. dipole->SetLocation(0., 0., -0.001);
  82. //dipole->SetLocation(-2.5,1.25,0);
  83. dipole->SetNumberOfFrequencies(1);
  84. dipole->SetFrequency(0, 2e7);
  85. //dipole->SetFrequency(1,2000);
  86. //dipole->SetPhase(0);
  87. //////////////////////////////////
  88. // H field Tests
  89. // // Define model, this is a difficult case
  90. // VectorXcr sigma(23);
  91. // sigma << 0., 0.0625, 0.0166666666666667, 0.04, 0.0181818181818182,
  92. // 0.05, 0.0222222222222222, 0.2, 0.025, 0.05, 0.025,
  93. // 0.05, 0.00833333333333333, 0.025, 0.05, 0.0222222222222222,
  94. // 0.0277777777777778, 0.0333333333333333, 0.0208333333333333,
  95. // 0.0181818181818182, 0.025, 0.04;
  96. //
  97. // VectorXr thick(21);
  98. // //thick << 10, 10, 10, 10, 10, 10;
  99. // thick << 11, 5, 5, 4, 3.5, 4.5, 4.5,
  100. // 5.5, 3, 36, 15, 3, 6, 3, 8,
  101. // 9, 10, 7, 8, 11, 68;
  102. //VectorXcr sigma(8);
  103. // sigma << 0., Complex(.01,0.01), 0.001, .005, .0034, .0023, .0024, .0001;
  104. VectorXcr sigma(8);
  105. sigma << 0., 3e-4, 5e-4, 1e-2, .5, 5e-6, .03, .04;
  106. //VectorXcr sigma(2);
  107. // sigma << 0., .0;
  108. //VectorXr susl(3);
  109. // susl << 1., 1., 1.5;//, 5e-6, 5e-6, 5e-6;
  110. //VectorXr sush(3);
  111. // sush << 1., 1.01;//, 5e-6, 5e-6, 5e-6;
  112. //VectorXr sustau(2);
  113. // sustau << 1., 1.0;//, 5e-6, 5e-6, 5e-6;
  114. //VectorXr susb(2);
  115. // susb << 0., .0;//, 5e-6, 5e-6, 5e-6;
  116. VectorXr thick(6);
  117. thick << 10 , 10 , 10 , 10 , 10 , 10;
  118. LayeredEarthEM *earth = LayeredEarthEM::New();
  119. earth->SetNumberOfLayers(8);
  120. //earth->SetNumberOfLayers(23);
  121. earth->SetLayerConductivity(sigma);
  122. //earth->SetLayerLowFreqSusceptibility(susl);
  123. //earth->SetLayerHighFreqSusceptibility(sush);
  124. //earth->SetLayerTauSusceptibility(sustau);
  125. //earth->SetLayerBreathSusceptibility(susb);
  126. //earth->SetLayerHighFreqSusceptibility(susl);
  127. //earth->SetLayerHighFreqSusceptibility(susl);
  128. earth->SetLayerThickness(thick);
  129. // just a test
  130. timer.end();
  131. // Receivers
  132. ReceiverPoints *receivers = ReceiverPoints::New();
  133. Vector3r loc;
  134. Real ox = 50.;
  135. Real oy = 70. ;
  136. Real depth = -40.13423;
  137. Real depth2 = depth;
  138. Real dz = .25;
  139. int nz = 3000;
  140. Real dx = .0;
  141. receivers->SetNumberOfReceivers(nz);
  142. int ir = 0;
  143. for (int iz=0; iz<nz; ++iz) {
  144. loc << ox, oy, depth;
  145. receivers->SetLocation(ir, loc);
  146. depth += dz;
  147. ox += dx;
  148. ++ ir;
  149. }
  150. // EmEarth
  151. EMEarth1D *EmEarth = EMEarth1D::New();
  152. EmEarth->AttachDipoleSource(dipole);
  153. EmEarth->AttachLayeredEarthEM(earth);
  154. EmEarth->AttachReceiverPoints(receivers);
  155. //EmEarth->SetFieldsToCalculate(BOTH);
  156. EmEarth->SetFieldsToCalculate(H);
  157. //EmEarth->SetHankelTransformMethod(GAUSSIANQUADRATURE);
  158. //EmEarth->SetHankelTransformMethod(DIGITALFILTERING);
  159. //EmEarth->SetHankelTransformMethod(FHTKEY201);
  160. //EmEarth->SetHankelTransformMethod(QWEKEY);
  161. EmEarth->SetHankelTransformMethod(FHTKEY51);
  162. timer.begin();
  163. EmEarth->MakeCalc3();
  164. //std::cout << "H\n" << receivers->GetHfield(0)
  165. // << "End of H" << std::endl;
  166. Real lemmaTime = timer.end();
  167. std::cout << lemmaTime << "\t";
  168. //std::cout << "Make Calc3\n";
  169. //receivers->ClearFields();
  170. //EmEarth->MakeCalc3();
  171. //std::cout << "Make Calc3 finished\n";
  172. #ifdef LEMMAUSEVTK
  173. int iif=0;
  174. vtkDataObject *_dataObjectxr = receivers->GetVtkDataObject(HFIELDREAL, iif,
  175. 0, receivers->GetNumberOfReceivers(), XCOMPONENT, ZCOORD);
  176. vtkDataObject *_dataObjectxi = receivers->GetVtkDataObject(HFIELDIMAG, iif,
  177. 0, receivers->GetNumberOfReceivers(), XCOMPONENT, ZCOORD);
  178. vtkDataObject *_dataObjectyr= receivers->GetVtkDataObject(HFIELDREAL, iif,
  179. 0, receivers->GetNumberOfReceivers(), YCOMPONENT, ZCOORD);
  180. vtkDataObject *_dataObjectyi = receivers->GetVtkDataObject(HFIELDIMAG, iif,
  181. 0, receivers->GetNumberOfReceivers(), YCOMPONENT, ZCOORD);
  182. vtkDataObject *_dataObjectzr = receivers->GetVtkDataObject(HFIELDREAL, iif,
  183. 0, receivers->GetNumberOfReceivers(), ZCOMPONENT, ZCOORD);
  184. vtkDataObject *_dataObjectzi = receivers->GetVtkDataObject(HFIELDIMAG, iif,
  185. 0, receivers->GetNumberOfReceivers(), ZCOMPONENT, ZCOORD);
  186. #endif
  187. /*
  188. // Write out to file
  189. depth = depth2;
  190. std::fstream real("reale_lay.dat", std::ios::out);
  191. std::fstream imag("image_lay.dat", std::ios::out);
  192. std::fstream hreal("real_lay.dat", std::ios::out);
  193. std::fstream himag("imag_lay.dat", std::ios::out);
  194. for (int iz=0; iz<nz; ++iz) {
  195. Vector3cr temp = receivers->GetEfield(0,iz);
  196. //std::cout << temp << endl;
  197. real << ox << "\t" << oy << "\t" << depth << "\t"
  198. << temp(0).real() << "\t" << temp(1).real()
  199. << "\t" << temp(2).real() << std::endl;
  200. imag << ox << "\t" << oy << "\t" << depth << "\t"
  201. << std::imag(temp(0)) << "\t" << std::imag(temp(1))
  202. << "\t" << std::imag(temp(2)) << std::endl;
  203. temp = receivers->GetHfield(0, iz);
  204. hreal << ox << "\t" << oy << "\t" << depth << "\t"
  205. << std::real(temp(0)) << "\t" << std::real(temp(1))
  206. << "\t" << std::real(temp(2)) << std::endl;
  207. himag << ox << "\t" << oy << "\t" << depth << "\t"
  208. << std::imag(temp(0)) << "\t" << std::imag(temp(1))
  209. << "\t" << std::imag(temp(2)) << std::endl;
  210. depth += dz;
  211. }
  212. real.close();
  213. imag.close();
  214. */
  215. // TODO check all field components for both E and H,
  216. // diplay in some sort of reasonable manner.
  217. // Then, and only then, optimise.
  218. receivers->ClearFields();
  219. //dipole->SetPolarisation(XPOLARISATION);
  220. //std::cout << "Fortran Calc" << std::endl;
  221. #ifdef KIHALEE_EM1D
  222. timer.begin();
  223. EmEarth->MakeCalc();
  224. Real fortranTime = timer.end();
  225. std::cout << fortranTime << "\t";
  226. std::cout << fortranTime/lemmaTime << std::endl;
  227. #endif
  228. #ifdef LEMMAUSEVTK
  229. vtkDataObject *_fdataObjectxr = receivers->GetVtkDataObject(HFIELDREAL, 0,
  230. 0, receivers->GetNumberOfReceivers(), XCOMPONENT, ZCOORD);
  231. vtkDataObject *_fdataObjectxi = receivers->GetVtkDataObject(HFIELDIMAG, 0,
  232. 0, receivers->GetNumberOfReceivers(), XCOMPONENT, ZCOORD);
  233. vtkDataObject *_fdataObjectyr = receivers->GetVtkDataObject(HFIELDREAL, 0,
  234. 0, receivers->GetNumberOfReceivers(), YCOMPONENT, ZCOORD);
  235. vtkDataObject *_fdataObjectyi = receivers->GetVtkDataObject(HFIELDIMAG, 0,
  236. 0, receivers->GetNumberOfReceivers(), YCOMPONENT, ZCOORD);
  237. vtkDataObject *_fdataObjectzr = receivers->GetVtkDataObject(HFIELDREAL, 0,
  238. 0, receivers->GetNumberOfReceivers(), ZCOMPONENT, ZCOORD);
  239. vtkDataObject *_fdataObjectzi = receivers->GetVtkDataObject(HFIELDIMAG, 0,
  240. 0, receivers->GetNumberOfReceivers(), ZCOMPONENT, ZCOORD);
  241. #endif
  242. // Simple 2D plot in vtk to compare algorithms
  243. #ifdef LEMMAUSEVTK
  244. vtkRenderer *_ren = vtkRenderer::New();
  245. vtkRenderWindow *_renwin = vtkRenderWindow::New();
  246. _renwin->AddRenderer(_ren);
  247. ///////////////////////////////////////////////
  248. // X Component Plot
  249. vtkXYPlotActor *_xplot = vtkXYPlotActor::New();
  250. //vtkXYPlotWidget *_wplot = vtkXYPlotWidget::New();
  251. // Add the datasets
  252. _xplot->AddDataObjectInput( _dataObjectxr);
  253. _xplot->AddDataObjectInput( _dataObjectxi);
  254. _xplot->AddDataObjectInput(_fdataObjectxr);
  255. _xplot->AddDataObjectInput(_fdataObjectxi);
  256. _xplot->SetTitle("H_x");
  257. _xplot->SetXTitle("");
  258. _xplot->SetYTitle("");
  259. _xplot->SetXValuesToValue();
  260. //_plot->SetXValuesToIndex();
  261. // set which parts of the data object are to be used for which axis
  262. _xplot->SetDataObjectXComponent(0,1);
  263. _xplot->SetDataObjectYComponent(0,0);
  264. _xplot->SetDataObjectXComponent(1,1);
  265. _xplot->SetDataObjectYComponent(1,0);
  266. _xplot->SetDataObjectXComponent(2,1);
  267. _xplot->SetDataObjectYComponent(2,0);
  268. _xplot->SetDataObjectXComponent(3,1);
  269. _xplot->SetDataObjectYComponent(3,0);
  270. _xplot->SetNumberOfXLabels(3);
  271. _xplot->SetPlotColor(0,1,1,1);
  272. _xplot->SetPlotColor(1,0,0,1);
  273. _xplot->SetPlotColor(2,1,0,0);
  274. _xplot->SetPlotColor(3,0,1,0);
  275. _xplot->SetPlotLabel(0, "Lemma real" );
  276. _xplot->SetPlotLabel(1, "Lemma imag" );
  277. _xplot->SetPlotLabel(2, "KiHa real" );
  278. _xplot->SetPlotLabel(3, "KiHa imag" );
  279. _xplot->GetProperty()->SetLineWidth(1);
  280. _xplot->GetProperty()->SetPointSize(22);
  281. _xplot->LegendOn();
  282. // _plot->SetPlotPoints(0,0);
  283. // _plot->SetPlotPoints(1,0);
  284. // _plot->SetPlotPoints(2,0);
  285. // _plot->SetPlotPoints(3,0);
  286. // _plot->SetPlotLines(0,0);
  287. // _plot->SetPlotLines(1,0);
  288. // _plot->SetPlotLines(2,0);
  289. // _plot->SetPlotLines(3,0);
  290. _xplot->PlotCurvePointsOff();
  291. _xplot->PlotCurveLinesOff();
  292. _xplot->GetPositionCoordinate()->SetValue(0.0, 0.67, 0);
  293. _xplot->GetPosition2Coordinate()->SetValue(1.0, 0.33, 0);
  294. //_plot->SetReverseYAxis(1); // Just flips axis, not data!
  295. ///////////////////////////////////////////////
  296. // Y Component Plot
  297. vtkXYPlotActor *_yplot = vtkXYPlotActor::New();
  298. //vtkXYPlotWidget *_wplot = vtkXYPlotWidget::New();
  299. // Add the datasets
  300. _yplot->AddDataObjectInput( _dataObjectyr);
  301. _yplot->AddDataObjectInput( _dataObjectyi);
  302. _yplot->AddDataObjectInput(_fdataObjectyr);
  303. _yplot->AddDataObjectInput(_fdataObjectyi);
  304. _yplot->SetTitle("H_y");
  305. _yplot->SetXTitle("");
  306. _yplot->SetYTitle("");
  307. _yplot->SetXValuesToValue();
  308. //_yplot->SetXValuesToIndex();
  309. // set which parts of the data object are to be used for which axis
  310. _yplot->SetDataObjectXComponent(0,1);
  311. _yplot->SetDataObjectYComponent(0,0);
  312. _yplot->SetDataObjectXComponent(1,1);
  313. _yplot->SetDataObjectYComponent(1,0);
  314. _yplot->SetDataObjectXComponent(2,1);
  315. _yplot->SetDataObjectYComponent(2,0);
  316. _yplot->SetDataObjectXComponent(3,1);
  317. _yplot->SetDataObjectYComponent(3,0);
  318. _yplot->SetNumberOfXLabels(3);
  319. _yplot->SetPlotColor(0,1,1,1);
  320. _yplot->SetPlotColor(1,0,0,1);
  321. _yplot->SetPlotColor(2,1,0,0);
  322. _yplot->SetPlotColor(3,0,1,0);
  323. _yplot->SetPlotLabel(0, "Lemma real" );
  324. _yplot->SetPlotLabel(1, "Lemma imag" );
  325. _yplot->SetPlotLabel(2, "KiHa real" );
  326. _yplot->SetPlotLabel(3, "KiHa imag" );
  327. _yplot->LegendOn();
  328. _yplot->GetProperty()->SetLineWidth(1);
  329. _yplot->GetProperty()->SetPointSize(12);
  330. _yplot->GetPositionCoordinate()->SetValue(0.00, 0.33, 0);
  331. _yplot->GetPosition2Coordinate()->SetValue(1.0, 0.33, 0);
  332. // _plot->SetPlotPoints(0,0);
  333. // _plot->SetPlotPoints(1,0);
  334. // _plot->SetPlotPoints(2,0);
  335. // _plot->SetPlotPoints(3,0);
  336. // _plot->SetPlotLines(0,0);
  337. // _plot->SetPlotLines(1,0);
  338. // _plot->SetPlotLines(2,0);
  339. // _plot->SetPlotLines(3,0);
  340. _yplot->PlotCurvePointsOff();
  341. _yplot->PlotCurveLinesOff();
  342. //_plot->SetReverseYAxis(1)
  343. ///////////////////////////////////////////////
  344. // Z Component Plot
  345. vtkXYPlotActor *_zplot = vtkXYPlotActor::New();
  346. //vtkXYPlotWidget *_wplot = vtkXYPlotWidget::New();
  347. // Add the datasets
  348. _zplot->AddDataObjectInput( _dataObjectzr);
  349. _zplot->AddDataObjectInput( _dataObjectzi);
  350. _zplot->AddDataObjectInput(_fdataObjectzr);
  351. _zplot->AddDataObjectInput(_fdataObjectzi);
  352. _zplot->SetTitle("H_z");
  353. _zplot->SetXTitle("Depth");
  354. _zplot->SetYTitle("");
  355. _zplot->SetXValuesToValue();
  356. //_plot->SetXValuesToIndex();
  357. // set which parts of the data object are to be used for which axis
  358. _zplot->SetDataObjectXComponent(0,1);
  359. _zplot->SetDataObjectYComponent(0,0);
  360. _zplot->SetDataObjectXComponent(1,1);
  361. _zplot->SetDataObjectYComponent(1,0);
  362. _zplot->SetDataObjectXComponent(2,1);
  363. _zplot->SetDataObjectYComponent(2,0);
  364. _zplot->SetDataObjectXComponent(3,1);
  365. _zplot->SetDataObjectYComponent(3,0);
  366. _zplot->SetNumberOfXLabels(3);
  367. _zplot->SetPlotColor(0,1,1,1);
  368. _zplot->SetPlotColor(1,0,0,1);
  369. _zplot->SetPlotColor(2,1,0,0);
  370. _zplot->SetPlotColor(3,0,1,0);
  371. _zplot->LegendOn();
  372. _zplot->SetPlotLabel(0, "Lemma real" );
  373. _zplot->SetPlotLabel(1, "Lemma imag" );
  374. _zplot->SetPlotLabel(2, "KiHa real" );
  375. _zplot->SetPlotLabel(3, "KiHa imag" );
  376. _zplot->GetProperty()->SetLineWidth(1);
  377. _zplot->GetProperty()->SetPointSize(14);
  378. _zplot->GetPositionCoordinate()->SetValue(0.0, 0.0, 0);
  379. _zplot->GetPosition2Coordinate()->SetValue(1.0, 0.33, 0);
  380. // _plot->SetPlotPoints(0,0);
  381. // _plot->SetPlotPoints(1,0);
  382. // _plot->SetPlotPoints(2,0);
  383. // _plot->SetPlotPoints(3,0);
  384. // _plot->SetPlotLines(0,0);
  385. // _plot->SetPlotLines(1,0);
  386. // _plot->SetPlotLines(2,0);
  387. // _plot->SetPlotLines(3,0);
  388. _zplot->PlotCurvePointsOn();
  389. _zplot->PlotCurveLinesOff();
  390. //_plot->SetReverseYAxis(1)
  391. _ren->AddActor2D(_xplot);
  392. _ren->AddActor2D(_yplot);
  393. _ren->AddActor2D(_zplot);
  394. vtkRenderWindowInteractor *_iren = vtkRenderWindowInteractor::New();
  395. _iren->SetRenderWindow(_renwin);
  396. _iren->Initialize();
  397. _renwin->Render();
  398. //_wplot->SetXYPlotActor(_plot);
  399. //_wplot->SetInteractor(_iren);
  400. //_wplot->EnabledOn();
  401. _iren->Start();
  402. // TODO clean up you dirty whore
  403. // _iren->Delete();
  404. // _plot->Delete();
  405. // _dataObject->Delete();
  406. // _fieldData->Delete();
  407. // _renwin->Delete();
  408. // _ren->Delete();
  409. // _xdatareal->Delete();
  410. // _ydata->Delete();
  411. // _depth->Delete();
  412. #endif
  413. EmEarth->Delete();
  414. receivers->Delete();
  415. earth->Delete();
  416. dipole->Delete();
  417. return 0;
  418. }