GPLIB++
angleavg.cpp
Go to the documentation of this file.
1 #include "VecMat.h"
2 #include "NetCDFTools.h"
3 #include <iostream>
4 #include <fstream>
5 #include <string>
6 #include <boost/algorithm/minmax_element.hpp>
7 #include <boost/algorithm/string/compare.hpp>
8 #include <boost/algorithm/string.hpp>
9 using namespace std;
10 using namespace gplib;
11 
12 int main()
13  {
14  string infilename, outfilename;
15 
16  string version = "$Id: 1dinvga.cpp 1136 2007-05-15 22:27:10Z max $";
17  cout << endl << endl;
18  cout << "Program " << version << endl;
19  cout
20  << "Create a two dimensional density plot from the models in the input file."
21  << endl;
22  cout << "This would usually be a concatenation of several '.plot' files "
23  << endl;
24  cout << "Specify the input file name and the output file name. " << endl;
25  cout << "Output is in netcdf format. " << endl;
26  cout << endl << endl;
27 
28  cout << "Input file: ";
29  cin >> infilename;
30 
31  ifstream infile(infilename.c_str());
32  std::vector<double> depth, anisotropy, angles;
33  double currdepth, curraniso, currangle;
34  while (infile.good())
35  {
36  infile >> currdepth >> curraniso >> currangle;
37  if (infile.good())
38  {
39  depth.push_back(currdepth);
40  anisotropy.push_back(curraniso);
41  angles.push_back(currangle);
42  }
43  }
44  if (depth.size() != parmvalue.size())
45  {
46  cerr << "Problem reading file !" << endl;
47  return 100;
48  }
49  else
50  {
51  cerr << "Read in file." << endl;
52  }
53  typedef vector<double>::const_iterator it;
54  pair<it, it> depthminmax =
55  boost::minmax_element(depth.begin(), depth.end());
56  pair<it, it> parmminmax = boost::minmax_element(parmvalue.begin(),
57  parmvalue.end());
58  const double depthstep = 1.0;
59  const double parmstep = 0.1;
60  const double deptheps = *(depthminmax.second) * 0.1;
61  const double parmeps = *(parmminmax.second) * 0.1;
62  const double mindepth = *(depthminmax.first);
63  const double maxdepth = *(depthminmax.second) + deptheps;
64  double minparm;
65  if (*(parmminmax.first) > 0.0)
66  minparm = *(parmminmax.first) * 0.9;
67  else
68  minparm = *(parmminmax.first) * 1.1;
69  const double maxparm = *(parmminmax.second) + parmeps;
70  const unsigned int ndepths = (maxdepth - mindepth) / depthstep;
71 
72  const unsigned int nparms = (maxparm - minparm) / parmstep + 1;
73  gplib::rmat Density(ndepths, nparms);
74 
75  fill_n(Density.data().begin(), ndepths * nparms, 0.0);
76  int nmodels = 0;
77  for (unsigned int i = 0; i < depth.size() - 1; ++i)
78  {
79  double nextdepth = depth.at(i + 1);
80  if (nextdepth < depth.at(i))
81  {
82  nextdepth = maxdepth;
83  }
84  unsigned int depthindex = (depth.at(i) - mindepth) / depthstep;
85  unsigned int nextdepthindex = (nextdepth - mindepth) / depthstep - 1;
86  unsigned int parmindex = (parmvalue.at(i) - minparm) / parmstep;
87  for (unsigned int j = depthindex; j < nextdepthindex; ++j)
88  Density(j, parmindex) += 1.0;
89  }
90 
91  gplib::rvec depthvals(ndepths);
92  gplib::rvec parmvals(nparms);
93  for (unsigned int i = 0; i < ndepths; ++i)
94  depthvals(i) = -1.0 * i * depthstep;
95  for (unsigned int i = 0; i < nparms; ++i)
96  {
97  if (wantlogparm)
98  {
99  parmvals(i) = log10(minparm) + (i + 0.5) * logparmstep;
100  }
101  else
102  {
103  parmvals(i) = minparm + (i + 0.5) * parmstep;
104  }
105  }
106  cout << "Output file: ";
107  cin >> outfilename;
108  ofstream avgfile((outfilename + ".avg").c_str());
109  for (unsigned int i = 0; i < ndepths; ++i) // for all depths
110  {
111  double avg = 0.0;
112  unsigned int count = 0;
113  for (unsigned int j = 0; j < nparms; ++j)
114  {
115  avg += Density(i, j) * parmvals(j);
116  count += Density(i, j);
117  }
118  avgfile << depthvals(i) << " " << avg / count << endl;
119  }
120  Write2DDataAsNetCDF(outfilename, Density, depthvals, parmvals, "Depth");
121 
122  }
123 
int main()
Definition: angleavg.cpp:12
string version
Definition: makeinput.cpp:16