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