Main Lemma Repository
Du kan inte välja fler än 25 ämnen Ämnen måste starta med en bokstav eller siffra, kan innehålla bindestreck ('-') och vara max 35 tecken långa.

EMDipEarth1D.cpp 3.9KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128
  1. // ===========================================================================
  2. //
  3. // Filename: utemdipearth1d.cpp
  4. //
  5. // Description:
  6. //
  7. // Version: 0.0
  8. // Created: 01/06/2010 04:27:56 PM
  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 <FDEM1D>
  34. using namespace Lemma;
  35. int main() {
  36. // Test with a single dipole
  37. auto dipole = DipoleSource::NewSP();
  38. dipole->SetType(GROUNDEDELECTRICDIPOLE);
  39. dipole->SetPolarisation(XPOLARISATION);
  40. dipole->SetNumberOfFrequencies(1);
  41. dipole->SetMoment(1);
  42. dipole->SetFrequency(0, 4400.1000);
  43. //dipole->SetPhase(0);
  44. //dipole->SetLocation( (VectorXr(3) << 49, 49, -1e-4).finished() );
  45. dipole->SetLocation( 49, 49, -1e-4 );
  46. // Define model
  47. VectorXcr sigma(2);
  48. sigma << 0., 1e-3;//, .1;//, .01, .001;
  49. VectorXr thick(1);
  50. thick << 10;//, 10, 10;
  51. auto earth = LayeredEarthEM::NewSP();
  52. earth->SetNumberOfLayers(2);
  53. earth->SetLayerConductivity(sigma);
  54. //earth->SetLayerThickness(thick);
  55. // Receivers
  56. auto receivers = FieldPoints::NewSP();
  57. Vector3r loc;
  58. //Real ox = -5.1456;
  59. //Real oy = 2.2350;
  60. Real ox = 50.;
  61. Real oy = 20.;
  62. Real depth = 18.10;
  63. Real dz = 2.6;
  64. int nz = 1;
  65. receivers->SetNumberOfPoints(nz);
  66. int ir = 0;
  67. for (int iz=0; iz<nz; ++iz) {
  68. loc << ox, oy, depth;
  69. //std::cout << "Receiver location " << loc.transpose() << std::endl;
  70. receivers->SetLocation(ir, loc);
  71. depth += dz;
  72. ++ ir;
  73. }
  74. //receivers->SetLocation(1, ox, oy, depth);
  75. auto EmEarth = EMEarth1D::NewSP();
  76. //EmEarth->SetHankelTransformMethod(DIGITALFILTERING);
  77. EmEarth->SetFieldsToCalculate(BOTH); // Fortran needs this
  78. EmEarth->AttachDipoleSource(dipole);
  79. EmEarth->AttachLayeredEarthEM(earth);
  80. EmEarth->AttachFieldPoints(receivers);
  81. //dipole->SetType(ELECTRICDIPOLE);
  82. // receivers->SetNumberOfReceivers(1);
  83. // //for
  84. // receivers->SetLocation(0, ox, oy, depth);
  85. std::cout << "Dipole location " << dipole->GetLocation( ).transpose() << "\n";
  86. std::cout << "Receiver location " << receivers->GetLocation(0).transpose() << "\n";
  87. std::cout << "C++\n";
  88. EmEarth->MakeCalc3();
  89. std::cout << receivers->GetEfield(0,0) << std::endl;
  90. //std::cout << receivers->GetEfield(0) << std::endl;
  91. //std::cout << receivers->GetHfield(1) << std::endl;
  92. //std::cout << receivers->GetEfield(1) << std::endl;
  93. //receivers->ClearFields();
  94. // swap tx rx posigion
  95. /*
  96. receivers->SetLocation(0, dipole->GetLocation());
  97. dipole->SetLocation(loc);
  98. EmEarth->MakeCalc3();
  99. std::cout << receivers->GetEfield(0,0) << std::endl;
  100. receivers->ClearFields();
  101. */
  102. auto lc = receivers->GetEfield(0,0);
  103. #ifdef KIHALEE_EM1D
  104. receivers->ClearFields();
  105. std::cout << "\nFORTRAN KiHa\n";
  106. EmEarth->MakeCalc();
  107. auto fc = receivers->GetEfield(0,0);
  108. // std::cout << receivers->GetHfield(0,0) << std::endl;
  109. std::cout << receivers->GetEfield(0,0) << std::endl;
  110. std::cout << "Difference norm |" << (lc - fc).norm() << "|" << std::endl;
  111. #endif
  112. }