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.

KernelV0.cpp 2.8KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980
  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() {
  21. auto earth = LayeredEarthEM::NewSP();
  22. earth->SetNumberOfLayers(3);
  23. earth->SetLayerConductivity( (VectorXcr(3) << Complex(0.,0), Complex(1./50.,0), Complex(1./100.)).finished() );
  24. earth->SetLayerThickness( (VectorXr(1) << 10).finished() );
  25. // Set mag field info
  26. // From NOAA, Laramie WY, June 9 2016, aligned with mag. north
  27. earth->SetMagneticFieldIncDecMag( 67, 0, 52750, NANOTESLA );
  28. // Transmitter loops
  29. auto Tx1 = CircularLoop(31, 15, 100, 100);
  30. auto Tx2 = CircularLoop(31, 15, 100, 120);
  31. //auto Tx1 = CircularLoop(60, 15, 0, 0); // was 60
  32. auto Kern = KernelV0::NewSP();
  33. Kern->PushCoil( "Coil 1", Tx1 );
  34. Kern->PushCoil( "Coil 2", Tx2 );
  35. Kern->SetLayeredEarthEM( earth );
  36. // std::cout << *Kern << std::endl;
  37. Kern->SetIntegrationSize( (Vector3r() << 200,200,2).finished() );
  38. Kern->SetIntegrationOrigin( (Vector3r() << 0,0,15).finished() );
  39. Kern->SetTolerance( 1e-13 );
  40. Kern->SetPulseDuration(0.020);
  41. Kern->SetPulseCurrent( VectorXr::LinSpaced( 20, .01, 200 ) ); // nbins, low, high
  42. Kern->SetDepthLayerInterfaces( VectorXr::LinSpaced( 20, .5, 50 ) );
  43. // We could, I suppose, take the earth model in here? For non-linear that
  44. // may be more natural to work with?
  45. std::vector<std::string> tx = {std::string("Coil 1")};
  46. std::vector<std::string> rx = {std::string("Coil 1")};
  47. Kern->CalculateK0( tx, rx, true );
  48. }
  49. std::shared_ptr<Lemma::PolygonalWireAntenna> CircularLoop ( int nd, Real Radius, Real Offsetx, Real Offsety ) {
  50. auto Tx1 = Lemma::PolygonalWireAntenna::NewSP();
  51. Tx1->SetNumberOfPoints(nd);
  52. VectorXr range = VectorXr::LinSpaced(nd, 0, 2*PI);
  53. int ii;
  54. for (ii=0; ii<nd; ++ii) {
  55. Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*std::cos(range(ii)), Offsety+Radius*std::sin(range(ii)), -1e-3));
  56. }
  57. //Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*1, Offsety, -1e-3));
  58. Tx1->SetCurrent(1.);
  59. Tx1->SetNumberOfTurns(1);
  60. Tx1->SetNumberOfFrequencies(1);
  61. Tx1->SetFrequency(0,2246);
  62. return Tx1;
  63. }