GPLIB++
MTStation.cpp
Go to the documentation of this file.
1 #include "MTStation.h"
2 #include "miscfunc.h"
3 //we only include the parsers for .j and .edi files
4 //if we have a proper installation of antlr
5 #ifdef HAVEANTLR
6 #include "EDILexer.hpp"
7 #include "EDIParser.hpp"
8 #include "JLexer.hpp"
9 #include "JParser.hpp"
10 #endif
11 //similarly we only activate reading netcdf files
12 //if netcdf is installed
13 #ifdef HAVENETCDF
14 #include "ReadWriteImpedances.h"
15 #endif
16 #include <iomanip>
17 #include <iterator>
18 #include <string>
19 #include "convert.h"
20 #include "FatalException.h"
21 #include <boost/filesystem/operations.hpp>
22 #include <boost/filesystem/convenience.hpp>
23 #include <boost/algorithm/string.hpp>
24 #include <boost/bind.hpp>
25 #include <boost/tokenizer.hpp>
26 #include <functional>
27 #include <fstream>
28 #include <cassert>
29 #include <map>
30 #include "Util.h"
31 using namespace std;
32 
33 namespace gplib
34  {
35  void TrimFilename(std::string &name)
36  {
37  boost::filesystem::path p(name);
38  name = p.stem().string();
39  }
40 
41  void MTStation::InitialSetup()
42  {
43  latitude = 0;
44  longitude = 0;
45  elevation = 0;
46  azimuth = 0;
47  }
48 
49  MTStation::MTStation()
50  {
51  InitialSetup();
52  }
53 
54  MTStation::~MTStation()
55  {
56  }
57 
58  MTStation::MTStation(const std::string filename)
59  {
60  InitialSetup();
61  GetData(filename);
62  }
63 
64  MTStation::MTStation(const int size) :
65  MTData(size)
66  {
67  InitialSetup();
68  Assign(size);
69  }
70  void MTStation::AssignAll(const int nfreq)
71  {
72  Assign(nfreq);
73  }
74  void MTStation::Assign(const int nfreq)
75  {
76  MTData.assign(nfreq, MTTensor());
77  TFData.assign(nfreq, MagneticTF());
78  }
79 
80  /*! The Update() method first updates all components individually and then
81  * recalculates all derived quantities that depend on more than one component
82  */
83 
84  void MTStation::Update()
85  {
86 
87  }
88 
89  /*! Rotates the component indicated by i by angle rotangle. Does not Update. */
90  void MTStation::DoRotate(const int i, const double angle)
91  {
92  dcomp newxx, newxy, newyx, newyy;
93  newxx = MTData.at(i).Zxx * std::pow(cos(angle), 2)
94  + (MTData.at(i).Zxy + MTData.at(i).Zyx) * sin(angle) * cos(angle)
95  + MTData.at(i).Zyy * std::pow(sin(angle), 2);
96  newxy = MTData.at(i).Zxy * std::pow(cos(angle), 2)
97  - (MTData.at(i).Zxx - MTData.at(i).Zyy) * sin(angle) * cos(angle)
98  - MTData.at(i).Zyx * std::pow(sin(angle), 2);
99  newyx = MTData.at(i).Zyx * std::pow(cos(angle), 2)
100  - (MTData.at(i).Zxx - MTData.at(i).Zyy) * sin(angle) * cos(angle)
101  - MTData.at(i).Zxy * std::pow(sin(angle), 2);
102  newyy = MTData.at(i).Zyy * std::pow(cos(angle), 2)
103  - (MTData.at(i).Zxy + MTData.at(i).Zyx) * sin(angle) * cos(angle)
104  + MTData.at(i).Zxx * std::pow(sin(angle), 2);
105  MTData.at(i).Zxx = newxx;
106  MTData.at(i).Zxy = newxy;
107  MTData.at(i).Zyx = newyx;
108  MTData.at(i).Zyy = newyy;
109  MTData.at(i).rotangle += angle;
110  }
111 
112  /*! Rotate(const double rotangle) rotates the impedance tensor by the angle rotangle in radian
113  * and updates all derived quantities
114  */
115  void MTStation::Rotate(const double rotangle)
116  {
117  for (unsigned int i = 0; i < MTData.size(); ++i)
118  {
119  DoRotate(i, rotangle);
120  }
121  Update();
122  }
123 
124  /*! Rotate() rotates the impedance tensor by the negative of the angle given in the rotangles field
125  */
126  void MTStation::Rotate()
127  {
128  for (unsigned int i = 0; i < MTData.size(); ++i)
129  {
130  DoRotate(i, -MTData.at(i).rotangle);
131  }
132  Update();
133  }
134  /*! Return a vector containing the available frequencies; This should be rewritten to be more efficient
135  * */
136  trealdata MTStation::GetFrequencies() const
137  {
138  trealdata temp(MTData.size());
139  transform(MTData.begin(), MTData.end(), temp.begin(), [](MTTensor val)
140  { return val.GetFrequency();});
141  return temp;
142  }
143 
144  void MTStation::SetFrequencies(const trealdata &freqs)
145  {
146  //assert(MTData.size() == freqs.size());
147  const unsigned int nfreqs = freqs.size();
148  if (MTData.size() != nfreqs) //this invalidates the data stored there before
149  {
150  MTData.resize(nfreqs);
151  TFData.resize(nfreqs);
152  }
153  for (unsigned int i = 0; i < nfreqs; ++i)
154  MTData.at(i).frequency = freqs.at(i);
155  assert(MTData.size() == TFData.size());
156 
157  }
158 
159  void MTStation::ReadMtf(const std::string filename)
160  {
161  std::ifstream mtffile(filename.c_str());
162  latitude = 0.0;
163  longitude = 0.0;
164  elevation = 0.0;
165  azimuth = 0.0;
166  name = filename;
167  std::string latline = FindToken(mtffile, ">LATITUDE");
168  if (mtffile.good())
169  {
170  typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
171  boost::char_separator<char> sep(":");
172  tokenizer tok(latline, sep);
173  tokenizer::iterator beg = tok.begin();
174  advance(beg, 1);
175  convert(*beg, latitude);
176  }
177  std::string lonline = FindToken(mtffile, ">LONGITUDE");
178  if (mtffile.good())
179  {
180  typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
181  boost::char_separator<char> sep(":");
182  tokenizer tok(lonline, sep);
183  tokenizer::iterator beg = tok.begin();
184  advance(beg, 1);
185  convert(*beg, longitude);
186  }
187  std::string line = FindToken(mtffile, "//SECTION=IMP");
188  double rotangle, period;
189  double zxx_re, zxx_im, zxy_re, zxy_im, zyx_re, zyx_im, zyy_re, zyy_im;
190  double dzxx, dzxy, dzyx, dzyy;
191 
192  while (mtffile.good())
193  {
194  MTTensor CurrentImp;
195  mtffile >> period >> rotangle >> zxx_re >> zxx_im >> dzxx >> zxy_re >> zxy_im
196  >> dzxy >> zyx_re >> zyx_im >> dzyx >> zyy_re >> zyy_im >> dzyy;
197  if (mtffile.good())
198  {
199 
200  CurrentImp.frequency = 1.0 / period;
201  CurrentImp.Zxx = (zxx_re - I * zxx_im);
202  CurrentImp.Zxy = (zxy_re - I * zxy_im);
203  CurrentImp.Zyx = (zyx_re - I * zyx_im);
204  CurrentImp.Zyy = (zyy_re - I * zyy_im);
205  CurrentImp.dZxx = dzxx;
206  CurrentImp.dZxy = dzxy;
207  CurrentImp.dZyx = dzyx;
208  CurrentImp.dZyy = dzyy;
209  CurrentImp.rotangle = rotangle;
210  MTData.push_back(CurrentImp);
211  TFData.push_back(MagneticTF());
212  }
213 
214  }
215 
216  }
217 
218  void MTStation::ReadZmm(const std::string filename)
219  {
220  std::ifstream zmmfile(filename.c_str());
221  latitude = 0.0;
222  longitude = 0.0;
223  elevation = 0.0;
224  azimuth = 0.0;
225  name = filename;
226  double dummy, rotangle;
227  std::string line = FindToken(zmmfile, "orientations");
228  zmmfile >> dummy >> rotangle;
229  try
230  {
231  while (zmmfile.good())
232  {
233  MTTensor CurrentImp;
234  double period = 0.0;
235  std::string periodline = FindToken(zmmfile, "period");
236  if (zmmfile.good())
237  {
238  typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
239 
240  boost::char_separator<char> sep(" ");
241  tokenizer tok(periodline, sep);
242  tokenizer::iterator beg = tok.begin();
243  advance(beg, 2);
244  std::cout << *beg << std::endl;
245  convert(*beg, period);
246  }
247  FindToken(zmmfile, "Transfer");
248  if (zmmfile.good())
249  {
250  double zxx_re, zxx_im, zxy_re, zxy_im, zyx_re, zyx_im, zyy_re, zyy_im;
251  zmmfile >> zxx_re >> zxx_im >> zxy_re >> zxy_im >> zyx_re >> zyx_im
252  >> zyy_re >> zyy_im;
253  CurrentImp.frequency = 1.0 / period;
254  CurrentImp.Zxx = (zxx_re + I * zxx_im);
255  CurrentImp.Zxy = (zxy_re + I * zxy_im);
256  CurrentImp.Zyx = (zyx_re + I * zyx_im);
257  CurrentImp.Zyy = (zyy_re + I * zyy_im);
258 
259  }
260  double s11_re, s11_im, s21_re, s21_im, s22_re, s22_im;
261  FindToken(zmmfile, "Inverse Coherent Signal Power Matrix");
262  if (zmmfile.good())
263  {
264  zmmfile >> s11_re >> s11_im >> s21_re >> s21_im >> s22_re >> s22_im;
265  }
266  FindToken(zmmfile, "Residual");
267  double cov11_re, cov11_im, cov21_re, cov21_im, cov22_re, cov22_im;
268  if (zmmfile.good())
269  {
270  zmmfile >> cov11_re >> cov11_im >> cov21_re >> cov21_im >> cov22_re
271  >> cov22_im;
272  }
273  CurrentImp.dZxx = sqrt(
274  std::abs(
275  complex<double>(cov11_re, cov11_im)
276  * complex<double>(s11_re, s11_im)));
277  CurrentImp.dZxy = sqrt(
278  std::abs(
279  complex<double>(cov11_re, cov11_im)
280  * complex<double>(s22_re, s22_im)));
281  CurrentImp.dZyx = sqrt(
282  std::abs(
283  complex<double>(cov22_re, cov22_im)
284  * complex<double>(s11_re, s11_im)));
285  CurrentImp.dZyy = sqrt(
286  std::abs(
287  complex<double>(cov22_re, cov22_im)
288  * complex<double>(s22_re, s22_im)));
289  CurrentImp.rotangle = rotangle;
290  MTData.push_back(CurrentImp);
291  TFData.push_back(MagneticTF());
292  }
293  } catch (...)
294  {
295 
296  }
297  }
298 
299 #ifdef HAVENETCDF
300  void MTStation::ReadNetCDF(const std::string filename)
301  {
302  std::vector<double> Frequencies,StatXCoord, StatYCoord, StatZCoord;
303  gplib::rvec Impedances;
304  ReadImpedancesFromNetCDF(filename, Frequencies, StatXCoord, StatYCoord, StatZCoord,Impedances);
305  if (StatXCoord.size() != 1)
306  {
307  throw FatalException("Cannot have data for more than one station in netcdf file !");
308  }
309  const size_t nfreq = Frequencies.size();
310  for (size_t i = 0; i < nfreq; ++i)
311  {
312  MTTensor CurrentImp;
313  CurrentImp.frequency = Frequencies[i];
314  CurrentImp.Zxx = (Impedances(i*8) + I * Impedances(i*8 +1));
315  CurrentImp.Zxy = (Impedances(i*8 +2) + I * Impedances(i*8 +3));
316  CurrentImp.Zyx = (Impedances(i*8 + 4) + I * Impedances(i*8 +5));
317  CurrentImp.Zyy = (Impedances(i*8 + 6) + I * Impedances(i*8 +7));
318  MTData.push_back(CurrentImp);
319  TFData.push_back(MagneticTF());
320  }
321  name = filename;
322  TrimFilename(name);
323  }
324 #endif
325 
326 #ifdef HAVEANTLR
327  void MTStation::ReadEdi(const std::string filename)
328  {
329  using namespace antlr;
330  std::ifstream infile(filename.c_str());
331  const double invalid = 1e30;
332 
333  if (infile)
334  {
335  try
336  {
337  EDILexer lexer(infile);
338  EDIParser parser(lexer);
339  parser.edi_file();
340  latitude = parser.latitude;
341  longitude = parser.longitude;
342  elevation = parser.elevation;
343  azimuth = parser.azimuth;
344  name = filename;
345  TrimFilename(name);
346  unsigned int valfreqs = 0;
347  for (unsigned int i = 0; i < parser.frequency.size(); ++i)
348  {
349  if (abs(parser.DataXX.at(i)) > invalid)
350  {
351  parser.frequency.at(i) = 0;
352  }
353  else
354  {
355  valfreqs++;
356  }
357  }
358  //frequency.assign(parser.frequency.size(),0);
359  //copy(parser.frequency.begin(),parser.frequency.end(),frequency.begin());
360  //rotangles.assign(parser.frequency.size(),0);
361  if (!parser.rotangles.empty())
362  {
363  azimuth += parser.rotangles.at(0) * 180. / PI;
364  }
365  Assign(valfreqs);
366  unsigned int mtindex = 0;
367  //There might be some invalid entries that we do not want to copy, we therefore
368  //go through all parsed entries with the for loop, but have a separate index
369  //for our storage class
370  for (unsigned int i = 0; i < parser.rotangles.size(); ++i)
371  {
372  if (abs(parser.DataXX.at(i)) < invalid)
373  {
374  MTData.at(mtindex).frequency = parser.frequency.at(i);
375  MTData.at(mtindex).rotangle = parser.rotangles.at(i)
376  * PI / 180.;
377  MTData.at(mtindex).Zxx = parser.DataXX.at(i);
378  MTData.at(mtindex).Zxy = parser.DataXY.at(i);
379  MTData.at(mtindex).Zyx = parser.DataYX.at(i);
380  MTData.at(mtindex).Zyy = parser.DataYY.at(i);
381  MTData.at(mtindex).dZxx = parser.dDataXX.at(i);
382  MTData.at(mtindex).dZxy = parser.dDataXY.at(i);
383  MTData.at(mtindex).dZyx = parser.dDataYX.at(i);
384  MTData.at(mtindex).dZyy = parser.dDataYY.at(i);
385  mtindex++;
386  }
387  }
388  mtindex = 0;
389  if (parser.DataZY.empty() == false)
390  {
391  for (unsigned int i = 0; i < parser.DataZY.size(); ++i)
392  {
393  if (abs(parser.DataZX.at(i)) < invalid)
394  {
395  TFData.at(mtindex).frequency = parser.frequency.at(
396  i);
397  TFData.at(mtindex).Tx = parser.DataZX.at(i);
398  TFData.at(mtindex).Ty = parser.DataZY.at(i);
399  TFData.at(mtindex).dTx = parser.dDataZX.at(i);
400  TFData.at(mtindex).dTy = parser.dDataZY.at(i);
401  mtindex++;
402  }
403  }
404  }
405  }
406  catch (ANTLRException& e)
407  {
408  cerr << "Parse exception: " << e.toString() << endl;
409  }
410  Update();
411  }
412  else
413  {
414  throw FatalException("File not found: " + filename);
415  }
416  }
417 #endif
418 
419  void MTStation::ReadPek1D(const std::string filename)
420  {
421  ifstream infile(filename.c_str());
422  if (infile)
423  {
424  const double convfactor = 1. / (1000. * mu);
425  double period, rxx, imxx, rxy, imxy, ryx, imyx, ryy, imyy;
426  while (infile.good())
427  {
428  infile >> period >> rxx >> imxx >> rxy >> imxy >> ryx >> imyx >> ryy
429  >> imyy;
430 
431  if (infile.good())
432  {
433  MTTensor CurrentImp;
434  CurrentImp.frequency = 1. / period;
435  CurrentImp.Zxx = (rxx + I * imxx) * convfactor;
436  CurrentImp.Zxy = (rxy + I * imxy) * convfactor;
437  CurrentImp.Zyx = (ryx + I * imyx) * convfactor;
438  CurrentImp.Zyy = (ryy + I * imyy) * convfactor;
439  MTData.push_back(CurrentImp);
440  TFData.push_back(MagneticTF());
441  }
442  }
443  name = filename;
444  TrimFilename(name);
445  }
446  else
447  {
448  throw FatalException("File not found: " + filename);
449  }
450  }
451 
452 #ifdef HAVEANTLR
453  void MTStation::ReadJ(const std::string filename)
454  {
455 
456  using namespace antlr;
457  ifstream infile(filename.c_str());
458 
459  if (infile)
460  {
461  try
462  {
463  JLexer lexer(infile);
464  JParser parser(lexer);
465  parser.jfile();
466  latitude = parser.latitude;
467  longitude = parser.longitude;
468  elevation = parser.elevation;
469  azimuth = parser.azimuth;
470  name = parser.name;
471  TrimFilename(name);
472  Assign(parser.frequency.size());
473 
474  //frequency = parser.frequency;
475  if (parser.zassigned == false && parser.rassigned == false)
476  {
477  cerr << "No MT data in file !" << endl;
478  }
479  else
480  {
481  for (unsigned int i = 0; i < parser.frequency.size(); ++i)
482  {
483  MTData.at(i).frequency = parser.frequency.at(i);
484  MTData.at(i).rotangle = parser.azimuth;
485  MTData.at(i).Zxx = parser.DataXX.at(i);
486  MTData.at(i).Zxy = parser.DataXY.at(i);
487  MTData.at(i).Zyx = parser.DataYX.at(i);
488  MTData.at(i).Zyy = parser.DataYY.at(i);
489  MTData.at(i).dZxx = parser.dDataXX.at(i);
490  MTData.at(i).dZxy = parser.dDataXY.at(i);
491  MTData.at(i).dZyx = parser.dDataYX.at(i);
492  MTData.at(i).dZyy = parser.dDataYY.at(i);
493  MTData.at(i).Rx = parser.Rx.at(i);
494  MTData.at(i).Ry = parser.Ry.at(i);
495  }
496  }
497  if (parser.tassigned)
498  {
499  for (unsigned int i = 0; i < parser.DataZY.size(); ++i)
500  {
501  TFData.at(i).frequency = parser.frequency.at(i);
502  TFData.at(i).Tx = parser.DataZX.at(i);
503  TFData.at(i).Ty = parser.DataZY.at(i);
504  TFData.at(i).dTx = parser.dDataZX.at(i);
505  TFData.at(i).dTy = parser.dDataZY.at(i);
506  TFData.at(i).Rz = parser.Rz.at(i);
507  }
508  }
509  }
510  catch (ANTLRException& e)
511  {
512  cerr << "Parse exception: " << e.toString() << endl;
513  }
514  Update();
515  }
516  else
517  {
518  throw FatalException("File not found: " + filename);
519  }
520 
521  }
522 #endif
523 
524  void MTStation::ReadMtt(const std::string filename)
525  {
526  ifstream infile;
527  double currentreal, currentimag;
528  int nentries = 0;
529  infile.open(filename.c_str());
530  while (infile.good())
531  {
532  infile >> currentreal;
533  ++nentries;
534  }
535  infile.close();
536  if (((nentries - 1) % 23) != 0)
537  throw FatalException("Number of records does not match expected: " + filename);
538  const int nrecords = (nentries - 1) / 23;
539  Assign(nrecords);
540  infile.open(filename.c_str());
541  int currentrecord = 0;
542  if (infile)
543  {
544  while (infile.good()) // read in inputfile
545  {
546  infile >> currentreal;
547  if (infile.good())
548  {
549  MTData.at(currentrecord).frequency = currentreal;
550  TFData.at(currentrecord).frequency = currentreal;
551  infile >> currentreal;
552  MTData.at(currentrecord).Nu = currentreal;
553  infile >> currentreal >> currentimag;
554  MTData.at(currentrecord).Zxx = (currentreal + I * currentimag);
555  infile >> currentreal >> currentimag;
556  MTData.at(currentrecord).Zxy = (currentreal + I * currentimag);
557  infile >> currentreal >> currentimag;
558  MTData.at(currentrecord).Zyx = (currentreal + I * currentimag);
559  infile >> currentreal >> currentimag;
560  MTData.at(currentrecord).Zyy = (currentreal + I * currentimag);
561  infile >> currentreal;
562  MTData.at(currentrecord).dZxx = currentreal;
563  infile >> currentreal;
564  MTData.at(currentrecord).dZxy = currentreal;
565  infile >> currentreal;
566  MTData.at(currentrecord).dZyx = currentreal;
567  infile >> currentreal;
568  MTData.at(currentrecord).dZyy = currentreal;
569  infile >> currentreal >> currentimag;
570  TFData.at(currentrecord).Tx = currentreal + I * currentimag;
571  infile >> currentreal >> currentimag;
572  TFData.at(currentrecord).Ty = currentreal + I * currentimag;
573  infile >> currentreal;
574  TFData.at(currentrecord).dTx = currentreal;
575  infile >> currentreal;
576  TFData.at(currentrecord).dTy = currentreal;
577  infile >> currentreal;
578  MTData.at(currentrecord).Rx = currentreal;
579  infile >> currentreal;
580  MTData.at(currentrecord).Ry = currentreal;
581  infile >> currentreal;
582  TFData.at(currentrecord).Rz = currentreal;
583  ++currentrecord;
584  }
585  }
586  infile.close();
587  name = filename;
588  TrimFilename(name);
589  }
590  else
591  {
592  throw FatalException("File not found: " + filename);
593  }
594  Update();
595  }
596 
597  void MTStation::GetData(const string filename)
598  {
599  if (!boost::filesystem::exists(filename))
600  throw FatalException("File does not exist : " + filename);
601  string ending = boost::filesystem::extension(filename);
602  boost::to_lower(ending);
603  if (ending == "")
604  {
605  cerr << "File has no extension ! " << filename << endl;
606  dataformat = unknown;
607  throw FatalException("File has no extension ! " + filename);
608  }
609  dataformat = MTFileTypes.at(ending);
610 
611  switch (dataformat)
612  {
613  case mtt:
614  ReadMtt(filename);
615  break;
616  case j:
617 #ifdef HAVEANTLR
618  ReadJ(filename);
619 #else
620  throw FatalException("Did not compile with .j support");
621 #endif
622  break;
623  case edi:
624 #ifdef HAVEANTLR
625  ReadEdi(filename);
626 #else
627  throw FatalException("Did not compile with .edi support");
628 #endif
629  break;
630  case pek:
631  ReadPek1D(filename);
632  break;
633  case zmm:
634  ReadZmm(filename);
635  break;
636  case mtf:
637  ReadMtf(filename);
638  break;
639  case nc:
640 #ifdef HAVENETCDF
641  ReadNetCDF(filename);
642 #else
643  throw FatalException("Did not compile with .nc support");
644 #endif
645  break;
646  default:
647  throw FatalException(
648  "Unknown File Format for Reading ! This should not happen");
649  }
650  assert(MTData.size() == TFData.size());
651  Rotate(); //rotate data to 0 degree
652 
653  }
654 
655  void MTStation::WriteData(const std::string filename)
656  {
657  name = filename;
658  WriteBack();
659  }
660 
661  void MTStation::WriteBack()
662  {
663  switch (dataformat)
664  {
665  case mtt:
666  WriteAsMtt(name.c_str());
667  break;
668  case j:
669  WriteAsJ(name.c_str());
670  break;
671  case edi:
672  WriteAsEdi(name.c_str());
673  break;
674  default:
675  throw FatalException(
676  "Unknown File Format for Writing ! This should not happen");
677  }
678 
679  }
680 
681  void MTStation::WriteAsMtt(const string filename)
682  {
683  Update();
684  WriteMtt((filename + ".mtt").c_str());
685  }
686 
687  void MTStation::WriteAsEdi(const string filename)
688  {
689  throw FatalException("WriteAsEdi not implemented yet !");
690  }
691 
692  void MTStation::WriteJBlock(boost::function<complex<double>(const MTTensor*)> Comp,
693  boost::function<double(const MTTensor*)> Err, ofstream &outfile,
694  const double convfactor)
695  {
696  for (unsigned int i = 0; i < MTData.size(); ++i)
697  {
698  if (MTData.at(i).frequency != 0)
699  {
700  outfile << setfill(' ') << setw(15) << resetiosflags(ios::fixed)
701  << 1. / MTData.at(i).frequency << " ";
702  outfile << setfill(' ') << setw(15)
703  << convfactor * boost::bind(Comp, &MTData.at(i))().real() << " ";
704  outfile << setfill(' ') << setw(15)
705  << convfactor * boost::bind(Comp, &MTData.at(i))().imag() << " ";
706  outfile << setfill(' ') << setw(15)
707  << convfactor * boost::bind(Err, &MTData.at(i))() << " ";
708  outfile << setfill(' ') << setw(15) << setiosflags(ios::fixed) << 1.000
709  << " " << endl;
710  }
711  }
712  }
713 
714  void MTStation::WriteAsJ(const string filename)
715  {
716  ofstream outfile;
717  outfile.open((filename + ".j").c_str());
718  int actfreqs = 0;
719  const double convfactor = 4. * acos(-1.) * 1e-4;
720 
721  outfile << "# J-File Produced by CJData.cpp" << endl;
722  outfile << ">LATITUDE = " << latitude << endl;
723  outfile << ">LONGITUDE = " << longitude << endl;
724  outfile << ">ELEVATION = " << elevation << endl;
725  outfile << ">AZIMUTH = " << azimuth << endl;
726  outfile << name << endl;
727  for (unsigned int i = 0; i < MTData.size(); ++i)
728  if (MTData.at(i).frequency != 0)
729  actfreqs++;
730  outfile << "ZXX S.I." << endl;
731  outfile << actfreqs << endl;
732  //WriteJBlock(DataXX,outfile,convfactor);
733  WriteJBlock(&MTTensor::GetZxx, &MTTensor::GetdZxx, outfile, convfactor);
734  outfile << "ZXY S.I." << endl;
735  outfile << actfreqs << endl;
736  WriteJBlock(&MTTensor::GetZxy, &MTTensor::GetdZxy, outfile, convfactor);
737 
738  outfile << "ZYX S.I." << endl;
739  outfile << actfreqs << endl;
740  WriteJBlock(&MTTensor::GetZyx, &MTTensor::GetdZyx, outfile, convfactor);
741 
742  outfile << "ZYY S.I." << endl;
743  outfile << actfreqs << endl;
744  WriteJBlock(&MTTensor::GetZyy, &MTTensor::GetdZyy, outfile, convfactor);
745 
746  /* if (norm(DataZY.Z.at(0)) !=0)
747  {
748  outfile << "TZX " << endl;
749  outfile << actfreqs << endl;
750  WriteJBlock(DataZX,outfile,1.);
751  }
752 
753  if (norm(DataZY.Z.at(0)) !=0)
754  {
755  outfile << "TZY " << endl;
756  outfile << actfreqs << endl;
757  WriteJBlock(DataZY,outfile,1.);
758  }*/
759  outfile.close();
760  }
761 
762  void inline MttLine(std::ofstream &outfile, double value)
763  {
764  outfile << setw(15) << setfill(' ') << setprecision(8) << setiosflags(ios::fixed)
765  << value << " ";
766  }
767 
768  void MTStation::WriteMtt(const std::string filename)
769  {
770  ofstream outfile;
771 
772  outfile.open(filename.c_str());
773  std::map<double, int> FreqOrder;
774  for (size_t i = 0; i < MTData.size(); ++i)
775  {
776  FreqOrder.insert(std::make_pair(MTData.at(i).frequency, i));
777  }
778  // index(frequency, mttindex);
779  assert(MTData.size() == TFData.size());
780  for (auto currfreq : FreqOrder) //write mtt-file
781  {
782  int i = currfreq.second;
783  if (MTData.at(i).frequency > 0)
784  {
785  //we write out the frequency in scientific notation
786  //to avoid loss of precision for low frequencies
787  outfile << setw(15) << setfill(' ') << setprecision(8)
788  << setiosflags(ios::scientific) << MTData.at(i).frequency;
789  outfile << " " << resetiosflags(ios::scientific) << setprecision(5)
790  << MTData.at(i).Nu << "\n";
791  //some of the older programs are very picky about
792  //how numbers are written, so we write the tensor
793  //elements in fixed format
794  MttLine(outfile, MTData.at(i).Zxx.real());
795  MttLine(outfile, MTData.at(i).Zxx.imag());
796  MttLine(outfile, MTData.at(i).Zxy.real());
797  MttLine(outfile, MTData.at(i).Zxy.imag());
798  MttLine(outfile, MTData.at(i).Zyx.real());
799  MttLine(outfile, MTData.at(i).Zyx.imag());
800  MttLine(outfile, MTData.at(i).Zyy.real());
801  MttLine(outfile, MTData.at(i).Zyy.imag());
802  outfile << "\n";
803  MttLine(outfile, MTData.at(i).dZxx);
804  MttLine(outfile, MTData.at(i).dZxy);
805  MttLine(outfile, MTData.at(i).dZyx);
806  MttLine(outfile, MTData.at(i).dZyy);
807  MttLine(outfile, TFData.at(i).Tx.real());
808  MttLine(outfile, TFData.at(i).Tx.imag());
809  MttLine(outfile, TFData.at(i).Ty.real());
810  MttLine(outfile, TFData.at(i).Ty.imag());
811  outfile << "\n";
812  MttLine(outfile, TFData.at(i).dTx);
813  MttLine(outfile, TFData.at(i).dTy);
814  MttLine(outfile, MTData.at(i).Rx);
815  MttLine(outfile, MTData.at(i).Ry);
816  MttLine(outfile, TFData.at(i).Rz);
817  outfile << "\n" << resetiosflags(ios::fixed);
818  //write out mtt file entries
819  }
820 
821  }
822  outfile.close();
823  }
824  MTStation::MTStation(const MTStation &old) :
825  latitude(old.latitude), longitude(old.longitude), elevation(old.elevation), name(
826  old.name), azimuth(old.azimuth), MTData(old.MTData), TFData(old.TFData),
827 
828  dataformat(old.dataformat)
829  {
830  Update();
831  }
832  MTStation& MTStation::operator=(const MTStation& source)
833  {
834  if (this == &source)
835  return *this;
836 
837  copy(source.MTData.begin(), source.MTData.end(), back_inserter(this->MTData));
838  copy(source.TFData.begin(), source.TFData.end(), back_inserter(this->TFData));
839 
840  this->latitude = source.latitude;
841  this->longitude = source.longitude;
842  this->elevation = source.elevation;
843  this->azimuth = source.azimuth;
844  this->name = source.name;
845  this->dataformat = source.dataformat;
846  this->Update();
847  return *this;
848  }
849  }
gplib::MTTensor::GetdZxx
double GetdZxx() const
Return tensor element errors.
Definition: MTTensor.h:152
gplib::MTStation::WriteData
virtual void WriteData(const std::string filename)
Definition: MTStation.cpp:655
gplib::MTStation::AssignAll
void AssignAll(const int nfreq)
Definition: MTStation.cpp:70
gplib::MTStation::operator=
MTStation & operator=(const MTStation &source)
Definition: MTStation.cpp:832
gplib::MTTensor::GetdZxy
double GetdZxy() const
Definition: MTTensor.h:156
gplib::MTStation::WriteAsEdi
void WriteAsEdi(const std::string filename)
Write data as edi (no functionality yet)
Definition: MTStation.cpp:687
gplib::MTStation::WriteAsJ
void WriteAsJ(const std::string filename)
Write data to j-file.
Definition: MTStation.cpp:714
gplib::MTTensor::GetZyx
std::complex< double > GetZyx() const
Definition: MTTensor.h:130
gplib::MTTensor::GetZxy
std::complex< double > GetZxy() const
Definition: MTTensor.h:126
gplib::TrimFilename
void TrimFilename(std::string &name)
Definition: MTStation.cpp:35
gplib::MttLine
void MttLine(std::ofstream &outfile, double value)
Definition: MTStation.cpp:762
gplib::MTStation
The class MTStation is used to store the transfer functions and related information for a MT-site...
Definition: MTStation.h:17
gplib::MTFileTypes
const std::map< std::string, MTStation::tmtdataformat > MTFileTypes
Definition: MTStation.h:157
MTData
CMTStation MTData
Definition: cadijoint.cpp:15
gplib::MagneticTF
Store th local magnetic transfer function (tipper)
Definition: MagneticTF.h:10
gplib::MTStation::GetData
virtual void GetData()
Definition: MTStation.h:141
gplib::MTTensor::GetdZyx
double GetdZyx() const
Definition: MTTensor.h:160
gplib::MTStation::Rotate
void Rotate(void)
Rotate to zero rotation angle.
Definition: MTStation.cpp:126
gplib::MTStation::WriteAsMtt
void WriteAsMtt(const std::string filename)
Write data in goettingen .mtt format.
Definition: MTStation.cpp:681
gplib::MTStation::Assign
void Assign(const int nfreq)
Assign() assigns zero to all derived quantities, this makes the calculation.
Definition: MTStation.cpp:74
gplib::MTTensor::GetdZyy
double GetdZyy() const
Definition: MTTensor.h:164
gplib::MTTensor
Stores MT-Tensor components at a single frequency, calculates derived quantities. ...
Definition: MTTensor.h:16
gplib::ReadImpedancesFromNetCDF
void ReadImpedancesFromNetCDF(const std::string &filename, std::vector< double > &Frequencies, std::vector< double > &StatXCoord, std::vector< double > &StatYCoord, std::vector< double > &StatZCoord, gplib::rvec &Impedances)
Read magnetotelluric impedances from a netcdf file.
Definition: ReadWriteImpedances.cpp:137
gplib::MTTensor::GetZxx
std::complex< double > GetZxx() const
Definition: MTTensor.h:122
size
const int size
Definition: perftest.cpp:14
gplib::MTStation::WriteBack
void WriteBack()
Write data back in original format, with filename given by station name.
Definition: MTStation.cpp:661
gplib::MTTensor::GetZyy
std::complex< double > GetZyy() const
Return tensor elements.
Definition: MTTensor.h:118
gplib::MTStation::GetFrequencies
trealdata GetFrequencies() const
return the available frequencies in a single vector
Definition: MTStation.cpp:136
gplib::MTStation::SetFrequencies
void SetFrequencies(const trealdata &freqs)
Set the frequencies of the tensor elements, invalidates the previously stored impedance data...
Definition: MTStation.cpp:144
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