Surface NMR forward modelling
您最多选择25个主题 主题必须以字母或数字开头,可以包含连字符 (-),并且长度不得超过35个字符

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124
  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 11/11/2016 02:44:37 PM
  11. * @version $Id$
  12. * @author Trevor Irons (ti)
  13. * @email tirons@egi.utah.edu
  14. * @copyright Copyright (c) 2016, University of Utah
  15. * @copyright Copyright (c) 2016, Lemma Software, LLC
  16. */
  17. #include <Merlin>
  18. using namespace Lemma;
  19. std::shared_ptr<PolygonalWireAntenna> CircularLoop ( int nd, Real radius, Real Offsetx, Real Offsety ) ;
  20. int main(int argc, char** argv) {
  21. if (argc < 6) { // 1 2 3 4 5
  22. std::cout << "./KernelV0 TxCoil.yaml RxCoil.yaml EMEarthModel.yaml AkvoDataSet.yaml Output.yaml \n";
  23. exit(EXIT_SUCCESS);
  24. }
  25. auto earth = LayeredEarthEM::NewSP();
  26. earth->SetNumberOfLayers(3);
  27. earth->SetLayerConductivity( (VectorXcr(3) << Complex(0.,0), Complex(1./50.,0), Complex(1./100.)).finished() );
  28. earth->SetLayerThickness( (VectorXr(1) << 10).finished() );
  29. // Set mag field info
  30. // From NOAA, Laramie WY, June 9 2016, aligned with mag. north
  31. earth->SetMagneticFieldIncDecMag( 67, 9, 52750, NANOTESLA );
  32. auto Tx1 = PolygonalWireAntenna::DeSerialize( YAML::LoadFile(argv[1]) );
  33. auto Tx2 = PolygonalWireAntenna::DeSerialize( YAML::LoadFile(argv[2]) );
  34. ////////////////////////////////////// _
  35. Tx1->SetCurrent(1.); // |
  36. Tx1->SetNumberOfTurns(1); // Set these from Akvo input!
  37. Tx1->SetNumberOfFrequencies(1); // |
  38. Tx1->SetFrequency(0,2246); // |
  39. Tx2->SetCurrent(1.); // |
  40. Tx2->SetNumberOfTurns(1); // \ | /
  41. Tx2->SetNumberOfFrequencies(1); // \ | /
  42. Tx2->SetFrequency(0,2246); // \ | /
  43. ////////////////////////////////////// _
  44. auto Kern = KernelV0::NewSP();
  45. Kern->PushCoil( "Coil 1", Tx1 );
  46. Kern->PushCoil( "Coil 2", Tx2 );
  47. Kern->SetLayeredEarthEM( earth );
  48. Kern->SetIntegrationSize( (Vector3r() << 200,200,200).finished() );
  49. Kern->SetIntegrationOrigin( (Vector3r() << 0,0,0).finished() );
  50. Kern->SetTolerance( 1e-12 ); // 1e-12
  51. auto AkvoDataNode = YAML::LoadFile(argv[4]);
  52. Kern->AlignWithAkvoDataset( AkvoDataNode );
  53. // These should to into AlignWithAkvoDataSet...
  54. Kern->SetPulseDuration( AkvoDataNode["pulseLength"][0].as<Real>() );
  55. Kern->SetPulseCurrent( AkvoDataNode["Pulses"]["Pulse 1"]["current"].as<VectorXr>() ); // nbins, low, high
  56. //VectorXr interfaces = VectorXr::LinSpaced( 41, .5, 45.5 ); // nlay, low, high
  57. VectorXr interfaces = VectorXr::LinSpaced( 51, .5, 45.5 ); // nlay, low, high
  58. Real thick = .5;
  59. for (int ilay=1; ilay<interfaces.size(); ++ilay) {
  60. interfaces(ilay) = interfaces(ilay-1) + thick;
  61. thick *= 1.05;
  62. }
  63. Kern->SetDepthLayerInterfaces( interfaces ); // nlay, low, high
  64. // We could, I suppose, take the earth model in here? For non-linear that
  65. // may be more natural to work with?
  66. //std::vector<std::string> tx = {std::string("Coil 1"), std::string("Coil 2") };
  67. std::vector<std::string> tx = {std::string("Coil 1")}; //, std::string("Coil 2") };
  68. std::vector<std::string> rx = {std::string("Coil 2")};
  69. Kern->CalculateK0( tx, rx, true );
  70. /*
  71. std::ofstream dout = std::ofstream(std::string("k-Tx2coil-Rx1coil-offset-")+ std::string(argv[1])+ std::string(".dat"));
  72. //std::ofstream dout = std::ofstream(std::string("k-coincident.dat"));
  73. dout << interfaces.transpose() << std::endl;
  74. dout << I.transpose() << std::endl;
  75. dout << "#real\n";
  76. dout << Kern->GetKernel().real() << std::endl;
  77. dout << "#imag\n";
  78. dout << Kern->GetKernel().imag() << std::endl;
  79. dout.close();
  80. */
  81. //std::ofstream out = std::ofstream(std::string("k-Tx2coil-Rx1coil-offset-")+std::string(argv[1])+std::string(".yaml"));
  82. std::ofstream out = std::ofstream(std::string(argv[5]));
  83. //std::ofstream out = std::ofstream(std::string("k-coincident.yaml"));
  84. out << *Kern;
  85. out.close();
  86. }
  87. std::shared_ptr<Lemma::PolygonalWireAntenna> CircularLoop ( int nd, Real Radius, Real Offsetx, Real Offsety ) {
  88. auto Tx1 = Lemma::PolygonalWireAntenna::NewSP();
  89. Tx1->SetNumberOfPoints(nd);
  90. VectorXr range = VectorXr::LinSpaced(nd, 0, 2*PI);
  91. int ii;
  92. for (ii=0; ii<nd; ++ii) {
  93. Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*std::cos(range(ii)), Offsety+Radius*std::sin(range(ii)), -1e-3));
  94. }
  95. //Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*1, Offsety, -1e-3));
  96. Tx1->SetCurrent(1.);
  97. Tx1->SetNumberOfTurns(1);
  98. Tx1->SetNumberOfFrequencies(1);
  99. Tx1->SetFrequency(0,2246);
  100. return Tx1;
  101. }