stackrf.cpp

Go to the documentation of this file.

#include "SeismicStationList.h"
#include "SeismicDataComp.h"
#include "VecMat.h"
#include "miscfunc.h"
#include <iostream>
#include <fstream>
#include <string>
#include "NumUtil.h"
#include "Util.h"

using namespace std;
using namespace gplib;

int main()
  {
    string version = "$Id: stackrf.cpp 1825 2009-11-03 13:38:06Z mmoorkamp $";
    cout << endl << endl;
    cout << "Program " << version << endl;
    cout << " Stack a number of receiver functions and";
    cout << " write out the averaged results." << endl;
    cout << " The program assumes that all RF have the same";
    cout << " sampling rate and shift, it does not correct for moveout."
        << endl;
    cout << endl << endl;
    string reclistname = AskFilename("File with list of filenames: ");
    string outfilename = AskFilename("Outfilename root: ");

    SeismicStationList RF;
    SeismicDataComp StackedRF;
    //read in the data from the filenames in reclistnames
    RF.ReadList(reclistname);
    //print some consistency information
    const unsigned int nrecfunc = RF.GetList().size();
    cout << "Read " << nrecfunc << " receiver functions." << endl;
    //the length of the first receiver function is the one we assume for everything
    const unsigned int rflength = RF.GetList().front()->GetData().size();
    //create objects for the average receiver function and the standard deviation
    gplib::rvec AvgRF(rflength), StdDev(rflength);
    fill_n(AvgRF.begin(), rflength, 0.0);
    fill_n(StdDev.begin(), rflength, 0.0);
    //go through all receiver functions
    for (unsigned int i = 0; i < nrecfunc; ++i)
      {
        //check that they have the same length
        if (RF.GetList().at(i)->GetData().size() != rflength)
          {
            cerr << "Incompatible data length in entry number " << i << endl;
            return 100;
          }
        transform(RF.GetList().at(i)->GetData().begin(),
            RF.GetList().at(i)->GetData().end(), AvgRF.begin(), AvgRF.begin(),
            std::plus<double>());
        //add current receiver function to average
        // for (unsigned int j = 0; j < rflength; ++j)
        // {
        //  AvgRF(j) += RF.GetList().at(i)->GetData().at(j);
        //  if (isnan(AvgRF(j)))
        //    cerr << i << " "<< j << " " << AvgRF(j) << std::endl;
        // }
      }
    //after everything has been summed, divide by the number of RF
    AvgRF *= 1. / nrecfunc;
    //and calculate the standard deviation
    for (unsigned int i = 0; i < nrecfunc; ++i)
      {
        for (unsigned int j = 0; j < rflength; ++j)
          {
            StdDev(j) += pow2(AvgRF(j) - RF.GetList().at(i)->GetData().at(
                j));
          }
      }
    //copy all header information from the first receiver functions
    StackedRF = *RF.GetList().front().get();
    //and copy the data from the averaged receiver function
    copy(AvgRF.begin(), AvgRF.end(), StackedRF.GetData().begin());
    //write out stacked RF in sac format
    StackedRF.WriteAsSac(outfilename + ".avg");
    StdDev *= 1. / (nrecfunc - 1);
    //write out stacked rf and standard deviation in ascii file
    ofstream outfile((outfilename + ".avg.asc").c_str());
    const double dt = RF.GetList().front()->GetDt();
    for (unsigned int i = 0; i < rflength; ++i)
      outfile << i * dt + RF.GetList().front()->GetB() << " " << AvgRF(i)
          << " " << sqrt(StdDev(i) / nrecfunc) << endl;
  }