params fixed for mkass
[Merlin.git] / examples / Coupling.cpp
1 /* This file is part of Lemma, a geophysical modelling and inversion API.
2 * More information is available at http://lemmasoftware.org
3 */
4
5 /* This Source Code Form is subject to the terms of the Mozilla Public
6 * License, v. 2.0. If a copy of the MPL was not distributed with this
7 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 */
9
10 /**
11 * @file
12 * @date 11/11/2016 02:44:37 PM
13 * @version $Id$
14 * @author Trevor Irons (ti)
15 * @email tirons@egi.utah.edu
16 * @copyright Copyright (c) 2016, University of Utah
17 * @copyright Copyright (c) 2016, Lemma Software, LLC
18 */
19
20 #include <Merlin>
21 using namespace Lemma;
22
23 static constexpr Real GAMMA = 2.67518e8; // MKS units
24
25 std::shared_ptr<PolygonalWireAntenna> CircularLoop ( int nd, Real radius, Real Offsetx, Real Offsety, Real wL ) ;
26 void MoveLoop( std::shared_ptr<PolygonalWireAntenna> Loop, int nd, Real Radius, Real Offsetx, Real Offsety, Real wL );
27
28 int main(int argc, char** argv) {
29 /*
30 if ( argc < 2 ) {
31 std::cout << "Calculates the coupling between two sNMR loops at the Larmor frequency. Usage\n"
32 << "\t./Coupling EarthModel.yaml" << std::endl;
33 exit(0);
34 }
35 //Real offset = atof(argv[1]);
36 auto earth = LayeredEarthEM::DeSerialize( YAML::LoadFile(argv[1]) );
37 */
38 // RedButtes model, also how you can generate your own files
39 auto earth = LayeredEarthEM::NewSP();
40 earth->SetNumberOfLayers(3);
41 //earth->SetLayerConductivity( (VectorXcr(3) << Complex(0.,0), Complex(1./50.,0), Complex(1./100.)).finished() );
42 earth->SetLayerConductivity( (VectorXcr(3) << Complex(0.,0), Complex(1./7.,0), Complex(1./100.)).finished() );
43 earth->SetLayerThickness( (VectorXr(1) << 10).finished() );
44 // Set mag field info
45 // From NOAA, Laramie WY, June 9 2016, aligned with mag. north
46 earth->SetMagneticFieldIncDecMag( 67, 0, 52750, NANOTESLA );
47 // auto sig = std::ofstream("SigmaModel.yaml");
48 // sig << *earth << std::endl;
49 // sig.close();
50
51 Real Larmor = earth->GetMagneticFieldMagnitude()*GAMMA/(2*PI);
52
53 // Transmitter loops
54 auto Tx1 = CircularLoop(21, 15, 50, 50, Larmor);
55 auto Tx2 = CircularLoop(21, 15, 50, 50, Larmor); // initially coincident
56
57 auto Kern = Coupling::NewSP();
58 Kern->PushCoil( "Coil 1", Tx1 );
59 Kern->PushCoil( "Coil 2", Tx2 );
60 Kern->SetLayeredEarthEM( earth );
61
62 Kern->SetIntegrationSize( (Vector3r() << 50,100,20).finished() );
63 Kern->SetIntegrationOrigin( (Vector3r() << 0,0,0.01).finished() );
64 Kern->SetTolerance( 1e-3 ); // 1e-12
65
66 std::vector<std::string> tx = {std::string("Coil 1")};//,std::string("Coil 2")};
67 std::vector<std::string> rx = {std::string("Coil 2")};
68 VectorXr Offsets = VectorXr::LinSpaced(6, 15.00, 23.0); // nbins, low, high
69
70 auto outfile = std::ofstream("coupling.dat");
71 for (int ii=0; ii< Offsets.size(); ++ii) {
72 MoveLoop(Tx2, 21, 15, 50, 50 + Offsets(ii), Larmor);
73 #ifdef LEMMAUSEVTK
74 Complex coupling = Kern->Calculate( tx, rx, true );
75 #else
76 Complex coupling = Kern->Calculate( tx, rx, false );
77 #endif
78 std::cout << "coupling " << coupling << std::endl;
79 outfile << Offsets(ii) << "\t" << std::real(coupling) << "\t" << std::imag(coupling) << std::endl;
80 }
81 outfile.close();
82 }
83
84 std::shared_ptr<Lemma::PolygonalWireAntenna> CircularLoop ( int nd, Real Radius, Real Offsetx, Real Offsety, Real wL ) {
85
86 auto Tx1 = Lemma::PolygonalWireAntenna::NewSP();
87 Tx1->SetNumberOfPoints(nd);
88
89 VectorXr range = VectorXr::LinSpaced(nd, 0, 2*PI);
90 int ii;
91 for (ii=0; ii<nd; ++ii) {
92 Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*std::cos(range(ii)), Offsety+Radius*std::sin(range(ii)), -1e-3));
93 }
94
95 Tx1->SetCurrent(1.);
96 Tx1->SetNumberOfTurns(1);
97 Tx1->SetNumberOfFrequencies(1);
98 Tx1->SetFrequency(0,wL);
99
100 return Tx1;
101 }
102
103 void MoveLoop( std::shared_ptr<Lemma::PolygonalWireAntenna> Tx1, int nd, Real Radius, Real Offsetx, Real Offsety, Real wL ) {
104
105 Tx1->SetNumberOfPoints(nd);
106
107 VectorXr range = VectorXr::LinSpaced(nd, 0, 2*PI);
108 int ii;
109 for (ii=0; ii<nd; ++ii) {
110 Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*std::cos(range(ii)), Offsety+Radius*std::sin(range(ii)), -1e-3));
111 }
112
113 Tx1->SetCurrent(1.);
114 Tx1->SetNumberOfTurns(1);
115 Tx1->SetNumberOfFrequencies(1);
116 Tx1->SetFrequency(0,wL);
117
118 }