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

Generated on Fri Jul 4 15:30:20 2008 for GPLIB++ by  doxygen 1.5.5