GPLIB++
C1DAnisoMTSynthData.cpp
Go to the documentation of this file.
1 #include "C1DAnisoMTSynthData.h"
2 #include <cassert>
3 #include <limits>
4 #include <fstream>
5 #include "NumUtil.h"
6 #include "FatalException.h"
7 
8 namespace gplib
9  {
10  //These periods are used as a default is no frequencies are specified
11  const int frequenzen = 50;
12  const float T[frequenzen] =
13  { 0.0025, 0.00308337, 0.00513901, 0.00625, 0.0104167, 0.0125, 0.0208333,
14  0.025, 0.0416667, 0.05, 0.0833333, 0.1, 0.166667, 0.2, 0.333333, 0.4,
15  0.667, 0.64, 1.06667, 1.29333, 2.15554, 2.59997, 4.33332, 5.19994,
16  8.66701, 10.666667, 16.000, 21.333, 32.000, 42.667, 64.0, 85.332,
17  127.992, 170.678, 256.016, 341.297, 512.033, 682.594, 1023.541,
18  1366.120, 2049.180, 2732.240, 4098.361, 5464.481, 8196.722,
19  10869.565, 16393.443, 28800, 43200, 86400 };
20 
21  void C1DAnisoMTSynthData::TransformRho()
22  {
23  assert(rho1.size() == rho2.size());
24  assert(rho1.size() == rho3.size());
25  assert(rho1.size() == strikes.size());
26  assert(rho1.size() == slants.size());
27  assert(rho1.size() == dips.size());
28 
29  const unsigned int nlayers = rho1.size();
30  effcond1.assign(nlayers, 0.0);
31  effcond2.assign(nlayers, 0.0);
32  effstrike.assign(nlayers, 0.0);
33  for (unsigned int i = 0; i < nlayers; ++i)
34  {
35  double s1 = 1. / rho1.at(i);//transform to conductivity
36  double s2 = 1. / rho2.at(i);
37  double s3 = 1. / rho3.at(i);
38  double strike = PI * strikes.at(i) / 180.0; //transform angles to radian
39  double dip = PI * dips.at(i) / 180.0;
40  double slant = PI * slants.at(i) / 180.0;
41  //calculate some frequently used sines and cosines
42  double sps = sin(strike);
43  double cps = cos(strike);
44  double sth = sin(dip);
45  double cth = cos(dip);
46  double sfi = sin(slant);
47  double cfi = cos(slant);
48  double pom1 = s1 * cfi * cfi + s2 * sfi * sfi;
49  double pom2 = s1 * sfi * sfi + s2 * cfi * cfi;
50  double pom3 = (s1 - s2) * sfi * cfi;
51  double c2ps = cps * cps;
52  double s2ps = sps * sps;
53  double c2th = cth * cth;
54  double s2th = sth * sth;
55  double csps = cps * sps;
56  double csth = cth * sth;
57  //contruct the conductivity tensor elements
58  double sg11 = pom1 * c2ps + pom2 * s2ps * c2th - 2. * pom3 * cth
59  * csps + s3 * s2th * s2ps;
60  double sg12 = pom1 * csps - pom2 * c2th * csps + pom3 * cth * (c2ps
61  - s2ps) - s3 * s2th * csps;
62  double sg13 = -pom2 * csth * sps + pom3 * sth * cps + s3 * csth
63  * sps;
64  double sg21 = sg12;
65  double sg22 = pom1 * s2ps + pom2 * c2ps * c2th + 2. * pom3 * cth
66  * csps + s3 * s2th * c2ps;
67  double sg23 = pom2 * csth * cps + pom3 * sth * sps - s3 * csth
68  * cps;
69  double sg31 = sg13;
70  double sg32 = sg23;
71  double sg33 = pom2 * s2th + s3 * c2th;
72  double axx = sg11 - sg13 * sg31 / sg33;
73  double axy = sg12 - sg13 * sg32 / sg33;
74  double ayx = sg21 - sg31 * sg23 / sg33;
75  double ayy = sg22 - sg23 * sg32 / sg33;
76  //calculate effective horizontal conductivity
77  double da12 = sqrt((axx - ayy) * (axx - ayy) + 4.0 * axy * ayx);
78  effcond1.at(i) = (0.5 * (axx + ayy + da12));
79  effcond2.at(i) = (0.5 * (axx + ayy - da12));
80  if (fcmp(da12, 0.0, std::numeric_limits<double>::epsilon()) != 0)
81  {
82  effstrike.at(i) = 0.5 * acos((axx - ayy) / da12);
83  }
84  else
85  {
86  effstrike.at(i) = 0.0;
87  }
88  if (axy < 0.0)
89  effstrike.at(i) *= -1.0;
90  //std::cout << effcond1.at(i) << " " << effcond2.at(i) << " " << effstrike.at(i)*180./PI << std::endl;
91 
92  }
93  }
94 
95  void C1DAnisoMTSynthData::GetData()
96  {
97 
98  calc_frequencies = GetFrequencies();
99  if (calc_frequencies.empty())
100  {
101  for (int i = 0; i < frequenzen; ++i)
102  calc_frequencies.push_back(1.0 / T[i]);
103  }
104  Assign(calc_frequencies.size());
105  TransformRho();
106  CalcZ();
107  Update();
108 
109  }
110 
111  inline dcomp dfp(const dcomp x)
112  {
113  //std::cout << 1.0 + exp(-2.0 * x) << std::endl;
114  return 1.0 + exp(-2.0 * x);
115  }
116 
117  inline dcomp dfm(const dcomp x)
118  {
119  //std::cout << 1.0 - exp(-2.0 * x) << std::endl;
120  return 1.0 - exp(-2.0 * x);
121  }
122 
123  void C1DAnisoMTSynthData::CalcZ()
124  {
125  dcomp omegamu;
126  const double convfactor = 1. / (1000. * mu);
127 
128  for (unsigned int i = 0; i < calc_frequencies.size(); ++i)
129  {
130  omegamu = -I * 8e-7 * PI * PI * calc_frequencies.at(i);
131  dcomp k0 = (1. - I) * 2.0e-3 * PI * sqrt(calc_frequencies.at(i)
132  / 10.0);
133  dcomp zxyrot = k0 / sqrt(effcond1.back());
134  dcomp zyxrot = -k0 / sqrt(effcond2.back());
135  double currstrike = effstrike.back();
136  MTTensor CurrZ(dcomp(0.0, 0.0), zxyrot, zyxrot, dcomp(0.0, 0.0),
137  calc_frequencies.at(i));
138  for (int layerindex = thicknesses.size() - 2; layerindex >= 0; --layerindex)
139  {
140  double currthick = 1000.0 * thicknesses.at(layerindex);
141  if (currstrike != effstrike.at(layerindex) && effcond1.at(layerindex)
142  != effcond2.at(layerindex))
143  {
144  CurrZ.Rotate(effstrike.at(layerindex) - currstrike);
145  currstrike = effstrike.at(layerindex);
146  }
147  MTTensor BottomZ(CurrZ);
148  //std::cout << CurrZ.GetZxy() * convfactor<< " " << CurrZ.GetZyx() * convfactor << std::endl;
149  dcomp dz1 = k0 / sqrt(effcond1.at(layerindex));
150  dcomp dz2 = k0 / sqrt(effcond2.at(layerindex));
151  dcomp ag1 = k0 * sqrt(effcond1.at(layerindex)) * currthick;
152  dcomp ag2 = k0 * sqrt(effcond2.at(layerindex)) * currthick;
153 
154  dcomp detzbottom = BottomZ.GetDet();
155  dcomp denominator = detzbottom * dfm(ag1) * dfm(ag2) / (dz1
156  * dz2) + BottomZ.GetZxy() * dfm(ag1) * dfp(ag2) / dz1
157  - BottomZ.GetZyx() * dfp(ag1) * dfm(ag2) / dz2 + dfp(ag1)
158  * dfp(ag2);
159  //std::cout << dz1 << dz2 << ag1 << ag2 << detzbottom << denominator << std::endl;
160  CurrZ.SetZxx() = 4.0 * BottomZ.GetZxx() * exp(-ag1 - ag2)
161  / denominator;
162  CurrZ.SetZxy() = (BottomZ.GetZxy() * dfp(ag1) * dfp(ag2)
163  - BottomZ.GetZyx() * dfm(ag1) * dfm(ag2) * dz1 / dz2
164  + detzbottom * dfp(ag1) * dfm(ag2) / dz2 + dfm(ag1) * dfp(
165  ag2) * dz1) / denominator;
166  CurrZ.SetZyx() = (BottomZ.GetZyx() * dfp(ag1) * dfp(ag2)
167  - BottomZ.GetZxy() * dfm(ag1) * dfm(ag2) * dz2 / dz1
168  - detzbottom * dfm(ag1) * dfp(ag2) / dz1 - dfp(ag1) * dfm(
169  ag2) * dz2) / denominator;
170  CurrZ.SetZyy() = 4.0 * BottomZ.GetZyy() * exp(-ag1 - ag2)
171  / denominator;
172  }
173  CurrZ.SetRotangle() = currstrike;
174  CurrZ.SetZxx() = convfactor * conj(CurrZ.SetZxx());
175  CurrZ.SetZxy() = convfactor * conj(CurrZ.SetZxy());
176  CurrZ.SetZyx() = convfactor * conj(CurrZ.SetZyx());
177  CurrZ.SetZyy() = convfactor * conj(CurrZ.SetZyy());
178  //std::cout << CurrZ.GetZxy() << " " << CurrZ.GetZyx() << std::endl << std::endl;
179  CurrZ.Rotate(-CurrZ.GetRotangle());
180  //
181  SetMTData().at(i) = CurrZ;
182  /*std::copy(effstrike.begin(),effstrike.end(),std::ostream_iterator<double>(std::cout," "));
183  std::cout << std::endl;
184  std::copy(effcond1.begin(),effcond1.end(),std::ostream_iterator<double>(std::cout," "));
185  std::cout << std::endl;
186  std::copy(effcond2.begin(),effcond2.end(),std::ostream_iterator<double>(std::cout," "));
187  std::cout << std::endl;*/
188 
189  }
190  }
191 
192  void C1DAnisoMTSynthData::WriteModel(std::string filename)
193  {
194  assert(rho1.size() == thicknesses.size());
195  assert(rho2.size() == thicknesses.size());
196  assert(rho3.size() == thicknesses.size());
197  assert(strikes.size() == thicknesses.size());
198  assert(slants.size() == thicknesses.size());
199  assert(dips.size() == thicknesses.size());
200  std::ofstream outfile(filename.c_str());
201  if (outfile)
202  {
203  outfile << thicknesses.size() << std::endl;
204  for (unsigned int i = 0; i < thicknesses.size(); ++i)
205  {
206  outfile << thicknesses.at(i) << " " << rho1.at(i) << " ";
207  outfile << rho2.at(i) << " " << rho3.at(i) << " ";
208  outfile << strikes.at(i) << " " << slants.at(i) << " ";
209  outfile << dips.at(i) << std::endl;
210  }
211  }
212  else
213  {
214  throw FatalException("Problem writing to file: " + filename);
215  }
216 
217  }
218 
219  void C1DAnisoMTSynthData::ReadModel(std::string filename)
220  {
221 
222  int layers;
223  double currthick, currho1, currho2, currho3, currstrike, currslant,
224  currdip;
225  int i = 0;
226 
227  std::ifstream infile(filename.c_str());
228  infile >> layers;
229  if (!infile.good() || layers < 0)
230  throw FatalException("Problem reading file: " + filename);
231  thicknesses.reserve(layers);
232  rho1.reserve(layers);
233  rho2.reserve(layers);
234  rho3.reserve(layers);
235  strikes.reserve(layers);
236  slants.reserve(layers);
237  dips.reserve(layers);
238  while (infile.good())
239  {
240  infile >> currthick >> currho1 >> currho2 >> currho3 >> currstrike
241  >> currslant >> currdip;
242  if (infile.good())
243  {
244  thicknesses.push_back(currthick);
245  rho1.push_back(currho1);
246  rho2.push_back(currho2);
247  rho3.push_back(currho3);
248  strikes.push_back(currstrike);
249  slants.push_back(currslant);
250  dips.push_back(currdip);
251  ++i;
252  }
253  }
254  assert(rho1.size() == thicknesses.size());
255  assert(rho2.size() == thicknesses.size());
256  assert(rho3.size() == thicknesses.size());
257  assert(strikes.size() == thicknesses.size());
258  assert(slants.size() == thicknesses.size());
259  assert(dips.size() == thicknesses.size());
260  assert(i==layers);
261  }
262 
263  void C1DAnisoMTSynthData::WritePlot(std::string filename)
264  {
265  std::ofstream outfile(filename.c_str());
266  double cumthick = 0;
267  TransformRho();
268  assert(effcond1.size() == thicknesses.size());
269  assert(effcond2.size() == thicknesses.size());
270  assert(effstrike.size() == thicknesses.size());
271  for (unsigned int i = 0; i < thicknesses.size(); ++i)
272  {
273  outfile << cumthick << " " << 1. / effcond1.at(i) << " ";
274  outfile << 1. / effcond2.at(i) << " " << effstrike.at(i) * 180.
275  / PI;
276  outfile << std::endl;
277  cumthick += thicknesses.at(i);
278  outfile << cumthick << " " << 1. / effcond1.at(i) << " ";
279  outfile << 1. / effcond2.at(i) << " " << effstrike.at(i) * 180.
280  / PI;
281  outfile << std::endl;
282  }
283  }
284 
285  gplib::rvec C1DAnisoMTSynthData::GetModelVector()
286  {
287  TransformRho();
288  assert(effcond1.size() == thicknesses.size());
289  assert(effcond2.size() == thicknesses.size());
290  assert(effstrike.size() == thicknesses.size());
291  const unsigned int length = thicknesses.size();
292  gplib::rvec result(length * 4);
293  for (unsigned int i = 0; i < length; ++i)
294  {
295  result(i) = log10(1. / effcond1.at(i));
296  result(i + length) = thicknesses.at(i);
297  result(i + 2 * length) = log10(effcond1.at(i) / effcond2.at(i));
298  result(i + 3 * length) = effstrike.at(i) * 180. / PI;
299  }
300  return result;
301  }
302 
303  C1DAnisoMTSynthData::C1DAnisoMTSynthData(const C1DAnisoMTSynthData &old) :
304  MTStation(old), calc_frequencies(old.calc_frequencies), thicknesses(
305  old.thicknesses), strikes(old.strikes), slants(old.slants), dips(
306  old.dips), rho1(old.rho1), rho2(old.rho2), rho3(old.rho3), effcond1(
307  old.effcond1), effcond2(old.effcond2), effstrike(old.effstrike)
308  {
309 
310  }
311  C1DAnisoMTSynthData::C1DAnisoMTSynthData()
312  {
313  }
314 
315  C1DAnisoMTSynthData::~C1DAnisoMTSynthData()
316  {
317  }
318  }
gplib::C1DAnisoMTSynthData::WritePlot
void WritePlot(std::string filename)
Definition: C1DAnisoMTSynthData.cpp:263
gplib::C1DAnisoMTSynthData::WriteModel
void WriteModel(std::string filename)
Definition: C1DAnisoMTSynthData.cpp:192
gplib::T
const float T[frequenzen]
Definition: C1DAnisoMTSynthData.cpp:12
gplib::C1DAnisoMTSynthData::GetData
virtual void GetData()
Calculate the synthetic data given the previously set parameters.
Definition: C1DAnisoMTSynthData.cpp:95
gplib::C1DAnisoMTSynthData
Calculate response of a 1D anisotropic model, code is based on Pek and Santos fortran code...
Definition: C1DAnisoMTSynthData.h:12
gplib::MTStation
The class MTStation is used to store the transfer functions and related information for a MT-site...
Definition: MTStation.h:17
gplib::dfp
dcomp dfp(const dcomp x)
Definition: C1DAnisoMTSynthData.cpp:111
gplib::dfm
dcomp dfm(const dcomp x)
Definition: C1DAnisoMTSynthData.cpp:117
gplib::frequenzen
const int frequenzen
Definition: C1DAnisoMTSynthData.cpp:11
gplib::C1DAnisoMTSynthData::ReadModel
void ReadModel(std::string filename)
Definition: C1DAnisoMTSynthData.cpp:219
gplib::MTStation::Assign
void Assign(const int nfreq)
Assign() assigns zero to all derived quantities, this makes the calculation.
Definition: MTStation.cpp:74
gplib::C1DAnisoMTSynthData::GetModelVector
gplib::rvec GetModelVector()
write model into file
Definition: C1DAnisoMTSynthData.cpp:285
gplib::MTStation::SetMTData
std::vector< MTTensor > & SetMTData()
Get the full vector of Tensor elements for reading and writing.
Definition: MTStation.h:124
gplib::MTStation::GetFrequencies
trealdata GetFrequencies() const
return the available frequencies in a single vector
Definition: MTStation.cpp:136
gplib::FatalException
The basic exception class for all errors that arise in gplib.
Definition: FatalException.h:11
gplib::MTStation::Update
void Update()
Update all derived quantities.
Definition: MTStation.cpp:84
Generated on Wed May 25 2016 14:21:06 for GPLIB++ by   doxygen 1.8.6