d52fe671426e834f81f9370c45c34d77b166354e
[Merlin.git] / examples / KernelV0-2.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
24 int main(int argc, char** argv) {
25
26 if (argc<4) {
27 std::cout << "./KernelV0-2 earth.yaml tx.yaml rx.yaml \n";
28 return(EXIT_SUCCESS);
29 }
30 std::cout << "Using earth model: " << argv[1] << std::endl;
31 auto earth = LayeredEarthEM::DeSerialize( YAML::LoadFile(argv[1]) );
32
33 std::cout << "Using transmitter: " << argv[2] << std::endl;
34 auto Tx = PolygonalWireAntenna::DeSerialize( YAML::LoadFile(argv[2]) );
35
36 std::cout << "Using receivers: " << argv[3] << std::endl;
37 auto Rx1 = PolygonalWireAntenna::DeSerialize( YAML::LoadFile(argv[3]) );
38
39 auto Kern = KernelV0::NewSP();
40 Kern->PushCoil( "Coil 1", Tx );
41 Kern->PushCoil( "Coil 2", Rx1 );
42 Kern->SetLayeredEarthEM( earth );
43
44 Kern->SetIntegrationSize( (Vector3r() << 20.2151538,20.438572,100).finished() );
45 Kern->SetIntegrationOrigin( (Vector3r() << -10, -10, .5).finished() );
46 Real tol(1e-13); // 13
47 Kern->SetTolerance( tol ); // 1e-12
48
49 // Kern->AlignWithAkvoDataset( YAML::LoadFile(argv[2]) );
50
51 Kern->SetPulseDuration(0.020);
52 VectorXr I(36);
53 // off from VC by 1.075926340216996
54 // Pulses from Wyoming Red Buttes exp 0
55 I << 397.4208916184016, 352.364477036168, 313.0112765842783, 278.37896394065376, 247.81424224324982,
56 220.77925043190442, 196.76493264105017, 175.31662279234038, 156.0044839325404, 138.73983004230124,
57 123.42064612625474, 109.82713394836259, 97.76534468972267, 87.06061858367781, 77.56000002944572, 69.1280780096311,
58 61.64250263640252, 54.99473044877554, 49.091182970515476, 43.84634004556388, 39.184136917167976, 35.03619319797924,
59 31.347205894128976, 28.06346770557137, 25.139117042424758, 22.53420773366429, 20.214205433283347,
60 18.144318026099942, 16.299965972298878, 14.652633628829891, 13.184271405688083, 11.870540177313893,
61 10.697057141915716, 9.64778948429609, 8.709338689612677, 7.871268012862094;
62 //Kern->SetPulseCurrent( VectorXr::LinSpaced( 1, 10, 200 ) ); // nbins, low, high
63 Kern->SetPulseCurrent( I ); // nbins, low, high
64
65 //VectorXr interfaces = VectorXr::LinSpaced( 41, .5, 45.5 ); // nlay, low, high
66 //VectorXr interfaces = VectorXr::LinSpaced( 61, .5, 45.5 ); // nlay, low, high
67 VectorXr interfaces = VectorXr::LinSpaced( 2, .5, 45.5 ); // nlay, low, high
68 Real thick = .1;
69 for (int ilay=1; ilay<interfaces.size(); ++ilay) {
70 interfaces(ilay) = interfaces(ilay-1) + thick;
71 thick *= 1.05;
72 }
73 Kern->SetDepthLayerInterfaces( interfaces ); // nlay, low, high
74
75 // We could, I suppose, take the earth model in here? For non-linear that
76 // may be more natural to work with?
77 std::vector<std::string> tx = {std::string("Coil 1")};
78 std::vector<std::string> rx = {std::string("Coil 2")};
79 Kern->CalculateK0( tx, rx, true ); // 3rd argument is vtk output
80
81 std::ofstream dout = std::ofstream(std::string("Rx-")+std::string(argv[3])+std::string(".dat"));
82 dout << "# Transmitters: ";
83 for (auto lp : tx) {
84 dout << lp << "\t";
85 }
86 dout << "\n# Receivers: ";
87 for (auto lp : rx) {
88 dout << lp << "\t";
89 }
90 dout << "\n# Tolerance: " << tol << std::endl;
91
92 dout << interfaces.transpose() << std::endl;
93 dout << I.transpose() << std::endl;
94 dout << "#real\n";
95 dout << Kern->GetKernel().real() << std::endl;
96 dout << "#imag\n";
97 dout << Kern->GetKernel().imag() << std::endl;
98 dout.close();
99
100 std::ofstream out = std::ofstream(std::string("Rx-")+std::string(argv[2])+std::string(".yaml"));
101 //std::ofstream out = std::ofstream(std::string("k-coincident.yaml"));
102 out << *Kern;
103 out.close();
104 }
105
106