angleavg.cpp
Go to the documentation of this file.00001 #include "VecMat.h"
00002 #include "NetCDFTools.h"
00003 #include <iostream>
00004 #include <fstream>
00005 #include <string>
00006 #include <boost/algorithm/minmax_element.hpp>
00007 #include <boost/algorithm/string/compare.hpp>
00008 #include <boost/algorithm/string.hpp>
00009 using namespace std;
00010 using namespace gplib;
00011
00012 int main()
00013 {
00014 string infilename, outfilename;
00015
00016 string version = "$Id: 1dinvga.cpp 1136 2007-05-15 22:27:10Z max $";
00017 cout << endl << endl;
00018 cout << "Program " << version << endl;
00019 cout
00020 << "Create a two dimensional density plot from the models in the input file."
00021 << endl;
00022 cout << "This would usually be a concatenation of several '.plot' files "
00023 << endl;
00024 cout << "Specify the input file name and the output file name. " << endl;
00025 cout << "Output is in netcdf format. " << endl;
00026 cout << endl << endl;
00027
00028 cout << "Input file: ";
00029 cin >> infilename;
00030
00031 ifstream infile(infilename.c_str());
00032 std::vector<double> depth, anisotropy, angles;
00033 double currdepth, curraniso, currangle;
00034 while (infile.good())
00035 {
00036 infile >> currdepth >> curraniso >> currangle;
00037 if (infile.good())
00038 {
00039 depth.push_back(currdepth);
00040 anisotropy.push_back(curraniso);
00041 angles.push_back(currangle);
00042 }
00043 }
00044 if (depth.size() != parmvalue.size())
00045 {
00046 cerr << "Problem reading file !" << endl;
00047 return 100;
00048 }
00049 else
00050 {
00051 cerr << "Read in file." << endl;
00052 }
00053 typedef vector<double>::const_iterator it;
00054 pair<it, it> depthminmax =
00055 boost::minmax_element(depth.begin(), depth.end());
00056 pair<it, it> parmminmax = boost::minmax_element(parmvalue.begin(),
00057 parmvalue.end());
00058 const double depthstep = 1.0;
00059 const double parmstep = 0.1;
00060 const double deptheps = *(depthminmax.second) * 0.1;
00061 const double parmeps = *(parmminmax.second) * 0.1;
00062 const double mindepth = *(depthminmax.first);
00063 const double maxdepth = *(depthminmax.second) + deptheps;
00064 double minparm;
00065 if (*(parmminmax.first) > 0.0)
00066 minparm = *(parmminmax.first) * 0.9;
00067 else
00068 minparm = *(parmminmax.first) * 1.1;
00069 const double maxparm = *(parmminmax.second) + parmeps;
00070 const unsigned int ndepths = (maxdepth - mindepth) / depthstep;
00071
00072 const unsigned int nparms = (maxparm - minparm) / parmstep + 1;
00073 gplib::rmat Density(ndepths, nparms);
00074
00075 fill_n(Density.data().begin(), ndepths * nparms, 0.0);
00076 int nmodels = 0;
00077 for (unsigned int i = 0; i < depth.size() - 1; ++i)
00078 {
00079 double nextdepth = depth.at(i + 1);
00080 if (nextdepth < depth.at(i))
00081 {
00082 nextdepth = maxdepth;
00083 }
00084 unsigned int depthindex = (depth.at(i) - mindepth) / depthstep;
00085 unsigned int nextdepthindex = (nextdepth - mindepth) / depthstep - 1;
00086 unsigned int parmindex = (parmvalue.at(i) - minparm) / parmstep;
00087 for (unsigned int j = depthindex; j < nextdepthindex; ++j)
00088 Density(j, parmindex) += 1.0;
00089 }
00090
00091 gplib::rvec depthvals(ndepths);
00092 gplib::rvec parmvals(nparms);
00093 for (unsigned int i = 0; i < ndepths; ++i)
00094 depthvals(i) = -1.0 * i * depthstep;
00095 for (unsigned int i = 0; i < nparms; ++i)
00096 {
00097 if (wantlogparm)
00098 {
00099 parmvals(i) = log10(minparm) + (i + 0.5) * logparmstep;
00100 }
00101 else
00102 {
00103 parmvals(i) = minparm + (i + 0.5) * parmstep;
00104 }
00105 }
00106 cout << "Output file: ";
00107 cin >> outfilename;
00108 ofstream avgfile((outfilename + ".avg").c_str());
00109 for (unsigned int i = 0; i < ndepths; ++i)
00110 {
00111 double avg = 0.0;
00112 unsigned int count = 0;
00113 for (unsigned int j = 0; j < nparms; ++j)
00114 {
00115 avg += Density(i, j) * parmvals(j);
00116 count += Density(i, j);
00117 }
00118 avgfile << depthvals(i) << " " << avg / count << endl;
00119 }
00120 Write2DDataAsNetCDF(outfilename, Density, depthvals, parmvals, "Depth");
00121
00122 }
00123