Ternary array modifications
[Merlin.git] / examples / KV0-3loops.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(int argc, char** argv) {
26
27 if (argc < 3) {
28 std::cout << "./KVo-3loops <offset> <tolerance> <rx>" << std::endl;
29 exit(0);
30 }
31
32 Real offset = atof(argv[1]);
33 std::cout << offset << std::endl;
34 Real tol = atof(argv[2]);
35
36 auto earth = LayeredEarthEM::NewSP();
37 earth->SetNumberOfLayers(3);
38 earth->SetLayerConductivity( (VectorXcr(3) << Complex(0.,0), Complex(1./50.,0), Complex(1./100.)).finished() );
39 earth->SetLayerThickness( (VectorXr(1) << 10).finished() );
40 // Set mag field info
41 // From NOAA, Laramie WY, June 9 2016, aligned with mag. north
42 earth->SetMagneticFieldIncDecMag( 67, 0, 52750, NANOTESLA );
43 //earth->SetMagneticFieldIncDecMag( 90, 0, 52750, NANOTESLA );
44 std::cout << "B0 " << earth->GetMagneticField( ).transpose() << std::endl;
45 std::cout << "hat BO " << earth->GetMagneticFieldUnitVector().transpose() << std::endl ;
46 std::cout << "hat |BO| " << earth->GetMagneticFieldUnitVector().norm() << std::endl ;
47
48 // Transmitter loops
49 auto Tx1 = CircularLoop(21, 15, 100+offset/2., 100-offset/2.);
50 auto Tx2 = CircularLoop(21, 15, 100+offset/2., 100+offset/2.);
51 auto Tx3 = CircularLoop(21, 15, 100-offset/2., 100 );
52
53
54
55 auto Kern = KernelV0::NewSP();
56 Kern->PushCoil( "Coil 1", Tx1 );
57 Kern->PushCoil( "Coil 2", Tx2 );
58 Kern->PushCoil( "Coil 3", Tx3 );
59 Kern->SetLayeredEarthEM( earth );
60
61 Kern->SetIntegrationSize( (Vector3r() << 200,200,200).finished() );
62 Kern->SetIntegrationOrigin( (Vector3r() << 0,0,0).finished() );
63 Kern->SetTolerance( tol ); // 1e-12
64
65 Kern->SetPulseDuration(0.020);
66 VectorXr I(36);
67
68 // off from VC by 1.075926340216996
69 // Pulses from Wyoming Red Buttes exp 0
70 I << 397.4208916184016, 352.364477036168, 313.0112765842783, 278.37896394065376, 247.81424224324982,
71 220.77925043190442, 196.76493264105017, 175.31662279234038, 156.0044839325404, 138.73983004230124,
72 123.42064612625474, 109.82713394836259, 97.76534468972267, 87.06061858367781, 77.56000002944572, 69.1280780096311,
73 61.64250263640252, 54.99473044877554, 49.091182970515476, 43.84634004556388, 39.184136917167976, 35.03619319797924,
74 31.347205894128976, 28.06346770557137, 25.139117042424758, 22.53420773366429, 20.214205433283347,
75 18.144318026099942, 16.299965972298878, 14.652633628829891, 13.184271405688083, 11.870540177313893,
76 10.697057141915716, 9.64778948429609, 8.709338689612677, 7.871268012862094;
77 //Kern->SetPulseCurrent( VectorXr::LinSpaced( 1, 10, 200 ) ); // nbins, low, high
78 Kern->SetPulseCurrent( I ); // nbins, low, high
79
80 //Kern->SetDepthLayerInterfaces( VectorXr::LinSpaced( 30, 3, 45.5 ) ); // nlay, low, high
81 VectorXr interfaces = VectorXr::LinSpaced( 41, .5, 45.5 ); // nlay, low, high
82 Real thick = .5;
83 for (int ilay=1; ilay<interfaces.size(); ++ilay) {
84 interfaces(ilay) = interfaces(ilay-1) + thick;
85 thick *= 1.05;
86 }
87 Kern->SetDepthLayerInterfaces( interfaces ); // nlay, low, high
88
89 // We could, I suppose, take the earth model in here? For non-linear that
90 // may be more natural to work with?
91 std::vector<std::string> tx = {std::string("Coil 1"), std::string("Coil 2"), std::string("Coil 3") };
92 std::vector<std::string> rx = {std::string("Coil 1"), std::string("Coil 3")};
93 //std::vector<std::string> rx = {std::string("Coil 1"), std::string("Coil 2"), std::string("Coil 3") };
94 //std::vector<std::string> rx = {std::string("Coil 1"), std::string("Coil 2")};//, std::string("Coil 3") };
95 //std::vector<std::string> rx = {std::string(argv[3])};
96 Kern->CalculateK0( tx, rx, false );
97
98 //std::ofstream dout = std::ofstream(std::string("k0-3Tx-RxCh-") + std::string(argv[3]) + std::string("-tol") + std::string(argv[1])+ std::string(".dat"));
99 std::ofstream dout = std::ofstream(std::string("k0-3Tx-RxCh-13") + std::string("-off-") + std::string(argv[1])+ std::string(".dat"));
100 dout << "# Transmitters: ";
101 for (auto lp : tx) {
102 dout << lp << "\t";
103 }
104 dout << "\n# Receivers: ";
105 for (auto lp : rx) {
106 dout << lp << "\t";
107 }
108 dout << "\n# Tolerance: " << tol << std::endl;
109 dout << "# Offset: " << offset << std::endl;
110 dout << "# Radius: " << 15 << std::endl;
111 dout << interfaces.transpose() << std::endl;
112 dout << Kern->GetPulseDuration()*I.transpose() << std::endl;
113 dout << "#real\n";
114 dout << Kern->GetKernel().real() << std::endl;
115 dout << "#imag\n";
116 dout << Kern->GetKernel().imag() << std::endl;
117 dout.close();
118
119 //std::ofstream out = std::ofstream(std::string("k0-3Tx-RxCh1-")+std::string(argv[1])+std::string(".yaml"));
120 std::ofstream out = std::ofstream(std::string("k0-3Tx-RxCh-13") + std::string("-off-") + std::string(argv[1])+ std::string(".yaml"));
121 //std::ofstream out = std::ofstream(std::string("k-coincident.yaml"));
122 out << *Kern;
123 out.close();
124 }
125
126 std::shared_ptr<Lemma::PolygonalWireAntenna> CircularLoop ( int nd, Real Radius, Real Offsetx, Real Offsety ) {
127
128 auto Tx1 = Lemma::PolygonalWireAntenna::NewSP();
129 Tx1->SetNumberOfPoints(nd);
130
131 VectorXr range = VectorXr::LinSpaced(nd, 0, 2*PI);
132 int ii;
133 for (ii=0; ii<nd; ++ii) {
134 Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*std::cos(range(ii)), Offsety+Radius*std::sin(range(ii)), -1e-3));
135 }
136 //Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*1, Offsety, -1e-3));
137
138 Tx1->SetCurrent(1.);
139 Tx1->SetNumberOfTurns(1);
140 Tx1->SetNumberOfFrequencies(1);
141 Tx1->SetFrequency(0,2246);
142
143 return Tx1;
144 }