Surface NMR forward modelling
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.

Coupling.cpp 4.3KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118
  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. static constexpr Real GAMMA = 2.67518e8; // MKS units
  20. std::shared_ptr<PolygonalWireAntenna> CircularLoop ( int nd, Real radius, Real Offsetx, Real Offsety, Real wL ) ;
  21. void MoveLoop( std::shared_ptr<PolygonalWireAntenna> Loop, int nd, Real Radius, Real Offsetx, Real Offsety, Real wL );
  22. int main(int argc, char** argv) {
  23. /*
  24. if ( argc < 2 ) {
  25. std::cout << "Calculates the coupling between two sNMR loops at the Larmor frequency. Usage\n"
  26. << "\t./Coupling EarthModel.yaml" << std::endl;
  27. exit(0);
  28. }
  29. //Real offset = atof(argv[1]);
  30. auto earth = LayeredEarthEM::DeSerialize( YAML::LoadFile(argv[1]) );
  31. */
  32. // RedButtes model, also how you can generate your own files
  33. auto earth = LayeredEarthEM::NewSP();
  34. earth->SetNumberOfLayers(3);
  35. //earth->SetLayerConductivity( (VectorXcr(3) << Complex(0.,0), Complex(1./50.,0), Complex(1./100.)).finished() );
  36. earth->SetLayerConductivity( (VectorXcr(3) << Complex(0.,0), Complex(1./7.,0), Complex(1./100.)).finished() );
  37. earth->SetLayerThickness( (VectorXr(1) << 10).finished() );
  38. // Set mag field info
  39. // From NOAA, Laramie WY, June 9 2016, aligned with mag. north
  40. earth->SetMagneticFieldIncDecMag( 67, 0, 52750, NANOTESLA );
  41. // auto sig = std::ofstream("SigmaModel.yaml");
  42. // sig << *earth << std::endl;
  43. // sig.close();
  44. Real Larmor = earth->GetMagneticFieldMagnitude()*GAMMA/(2*PI);
  45. // Transmitter loops
  46. auto Tx1 = CircularLoop(21, 15, 50, 50, Larmor);
  47. auto Tx2 = CircularLoop(21, 15, 50, 50, Larmor); // initially coincident
  48. auto Kern = Coupling::NewSP();
  49. Kern->PushCoil( "Coil 1", Tx1 );
  50. Kern->PushCoil( "Coil 2", Tx2 );
  51. Kern->SetLayeredEarthEM( earth );
  52. Kern->SetIntegrationSize( (Vector3r() << 50,100,20).finished() );
  53. Kern->SetIntegrationOrigin( (Vector3r() << 0,0,0.01).finished() );
  54. Kern->SetTolerance( 1e-12 ); // 1e-12
  55. std::vector<std::string> tx = {std::string("Coil 1")};//,std::string("Coil 2")};
  56. std::vector<std::string> rx = {std::string("Coil 2")};
  57. VectorXr Offsets = VectorXr::LinSpaced(6, 15.00, 23.0); // nbins, low, high
  58. auto outfile = std::ofstream("coupling.dat");
  59. for (int ii=0; ii< Offsets.size(); ++ii) {
  60. MoveLoop(Tx2, 21, 15, 50, 50 + Offsets(ii), Larmor);
  61. #ifdef LEMMAUSEVTK
  62. Complex coupling = Kern->Calculate( tx, rx, true );
  63. #else
  64. Complex coupling = Kern->Calculate( tx, rx, false );
  65. #endif
  66. std::cout << "coupling " << coupling << std::endl;
  67. outfile << Offsets(ii) << "\t" << std::real(coupling) << "\t" << std::imag(coupling) << std::endl;
  68. }
  69. outfile.close();
  70. }
  71. std::shared_ptr<Lemma::PolygonalWireAntenna> CircularLoop ( int nd, Real Radius, Real Offsetx, Real Offsety, Real wL ) {
  72. auto Tx1 = Lemma::PolygonalWireAntenna::NewSP();
  73. Tx1->SetNumberOfPoints(nd);
  74. VectorXr range = VectorXr::LinSpaced(nd, 0, 2*PI);
  75. int ii;
  76. for (ii=0; ii<nd; ++ii) {
  77. Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*std::cos(range(ii)), Offsety+Radius*std::sin(range(ii)), -1e-3));
  78. }
  79. Tx1->SetCurrent(1.);
  80. Tx1->SetNumberOfTurns(1);
  81. Tx1->SetNumberOfFrequencies(1);
  82. Tx1->SetFrequency(0,wL);
  83. return Tx1;
  84. }
  85. void MoveLoop( std::shared_ptr<Lemma::PolygonalWireAntenna> Tx1, int nd, Real Radius, Real Offsetx, Real Offsety, Real wL ) {
  86. Tx1->SetNumberOfPoints(nd);
  87. VectorXr range = VectorXr::LinSpaced(nd, 0, 2*PI);
  88. int ii;
  89. for (ii=0; ii<nd; ++ii) {
  90. Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*std::cos(range(ii)), Offsety+Radius*std::sin(range(ii)), -1e-3));
  91. }
  92. Tx1->SetCurrent(1.);
  93. Tx1->SetNumberOfTurns(1);
  94. Tx1->SetNumberOfFrequencies(1);
  95. Tx1->SetFrequency(0,wL);
  96. }