GPLIB++
rotts.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <string>
3 #include <numeric>
4 #include "TimeSeriesData.h"
5 #include "FatalException.h"
6 #include "statutils.h"
7 #include "Util.h"
8 
9 using namespace std;
10 using namespace gplib;
11 
12 string version = "$Id: rotts.cpp 1844 2010-04-12 11:34:25Z mmoorkamp $";
13 
14 /*!
15  * \addtogroup UtilProgs Utility Programs
16  *@{
17  * \file
18  * Rotate a MT time-series so that the By component has a zero mean.
19  */
20 
21 int main(int argc, char *argv[])
22  {
23  string infilename;
24 
25  cout << "This is rotts: Rotate time series to a 0 mean By component"
26  << endl << endl;
27  cout << " Usage: rotts filename1" << endl;
28  cout
29  << " The output files will have the same name as the input files + .rot "
30  << endl << endl;
31  cout << " This is Version: " << version << endl << endl;
32  if (argc == 2)
33  {
34  infilename = argv[1];
35  }
36  else
37  {
38  infilename = AskFilename("Input Filename: ");
39 
40  }
41  cout << "Rotating ... " << endl;
42  TimeSeriesData Data1;
43  try
44  {
45  Data1.GetData(infilename);
46 
47  double hymedian = Median(Data1.GetData().GetHy().GetData().begin(),
48  Data1.GetData().GetHy().GetData().end());
49  double hxmedian = Median(Data1.GetData().GetHx().GetData().begin(),
50  Data1.GetData().GetHx().GetData().end());
51  cout << "Hy Median before: " << hymedian << endl;
52 
53  double phi = hymedian / sqrt(std::pow(hxmedian, 2) + std::pow(hymedian, 2));
54  cout << "Estimated angle: " << phi << endl;
55  double newhx, newhy;
56  double sinphi = sin(-phi);
57  double cosphi = cos(-phi);
58  const unsigned int ndata = Data1.GetData().GetHy().GetData().size();
59  for (unsigned int i = 0; i < ndata; ++i)
60  {
61  newhx = cosphi * Data1.GetData().GetHx().GetData().at(i) + sinphi
62  * Data1.GetData().GetHy().GetData().at(i);
63  newhy = sinphi * Data1.GetData().GetHx().GetData().at(i) + cosphi
64  * Data1.GetData().GetHy().GetData().at(i);
65  Data1.GetData().GetHx().GetData().at(i) = newhx;
66  Data1.GetData().GetHy().GetData().at(i) = newhy;
67  }
68  hymedian = Median(Data1.GetData().GetHy().GetData().begin(),
69  Data1.GetData().GetHy().GetData().end());
70  cout << "Hy Median after: " << hymedian << endl;
71  Data1.WriteAsBirrp(infilename + ".rot");
72  } catch (const FatalException &e)
73  {
74  cerr << e.what() << endl;
75  }
76  cout << "... done " << endl;
77  }
78 /*@}*/
79 
std::iterator_traits< InputIterator >::value_type Median(InputIterator begin, InputIterator end)
Calculate the median for a vector style container.
Definition: statutils.h:102
std::vector< double > & GetData()
Access for data vector, for ease of use and efficiency we return a reference.
TimeSeries & GetData()
return a reference to the actual object stored in the pointer
string version
Definition: rotts.cpp:12
TimeSeriesComponent & GetHy()
Definition: TimeSeries.h:39
int main()
Definition: angleavg.cpp:12
void WriteAsBirrp(std::string filename_base)
Write data to file in ascii format for birrp processing.
TimeSeriesData stores a pointer to the different components of magnetotelluric data and provides func...
The basic exception class for all errors that arise in gplib.
TimeSeriesComponent & GetHx()
Access function for Hx, returns reference for efficiency.
Definition: TimeSeries.h:35