calcrec.cpp

Go to the documentation of this file.

#include <iostream>
#include <string>
#include <algorithm>
#include "Util.h"
#include "RecCalc.h"
#include "CalcRecConf.h"
#include "SeismicDataComp.h"
#include "FatalException.h"
#include "FilterFunc.h"
#include "statutils.h"

/*! \file
 * This program calculates a receiver function from single radial and vertical components.
 * When started the program only asks for the name of the files for these two components.
 * The output receiver function has the same name as the radial component with .rec appended.
 * All other aspects of the calculation are specified in the file calcrec.conf in the directory
 * in which the program is executed.
 *
 * The parameters for that file are as follows:
 *
 * <TABLE>
 * <TR> <TD> <B>Name </B></TD> <TD> <B>Description </B></TD></TR>
 * <TR> <TD> recmethod </TD>  <TD> Can be specdiv for spectral division method or iterdecon for iterative deconvolution. </TD></TR>
 * <TR> <TD> omega </TD>  <TD>The center of the gaussian used to filter the receiver function during estimation. </TD> </TR>
 * <TR> <TD> sigma </TD>  <TD>The width of the gaussian used to filter the receiver function during estimation.
 *    Larger values mean a more broad band reciever function that contains more detailed structure but is also susceptible to noise.</TD> </TR>
 * <TR> <TD> shift </TD>  <TD>The shift in seconds of the initial correlation peak. This is purely used for display purposes. </TD> </TR>
 * <TR> <TD> cc </TD>  <TD>(only specdiv) The waterlevel relative to the maximum autopower of the vertical component. </TD> </TR>
 * <TR> <TD> upperfreq </TD>  <TD>The upper frequency in Hz for the bandpass filter, if negative or not set it is ignored </TD> </TR>
 * <TR> <TD> lowerfreq </TD>  <TD>The lower frequency in Hz for the bandpass filter, if negative or not set it is ignored </TD> </TR>
 * </TABLE>
 * More details about these parameters can be found for example in
 * Langston, 1979, JGR 84(B9) 4749-4762
 * Liggoria, 1999, BSSA 89(5), 1395-1400
 */
using namespace std;
using namespace gplib;

string version = "$Id: calcrec.cpp 1839 2010-03-05 17:04:34Z mmoorkamp $";

int main(int argc, char *argv[])
  {
    // write some info
    cout << " This is datarec: Calculate receiver functions from input data"
        << endl;
    cout << " Reads in radial and vertical components in SAC format" << endl;
    cout
        << " Outputs a file with the name of the radial component and .rec appended"
        << endl;
    cout << " Some behaviour can be configured with the file calcrec.conf"
        << endl;
    cout << " This is Version: " << version << endl << endl;
    string modelfilename;

    CalcRecConf Config;
    SeismicDataComp Radial, Vertical;
    ResPkModel Model;
    string outfilename;
    try
      {
        //get the settings from the configuration file
        Config.GetData("calcrec.conf");
        //the default method is spectral division
        RecCalc::trfmethod rfmethod = RecCalc::specdiv;
        if (Config.recmethod == "iterdecon")
          {
            rfmethod = RecCalc::iterdecon;
          }
        string radfilename, verfilename;
        //get the names for the radial and vertical components
        //either from the command line or by asking the user
        if (argc == 3)
          {
            radfilename = argv[1];
            verfilename = argv[2];
          }
        else
          {
            radfilename = AskFilename("Radial Component: ");
            verfilename = AskFilename("Vertical Component: ");
          }
        //the receiver function name is determined from the radial component
        outfilename = radfilename + ".rec";
        //read in the two components
        Radial.ReadData(radfilename);
        Vertical.ReadData(verfilename);
        if (Config.upperfreq > 0.0 && Config.lowerfreq > 0.0
            && Config.upperfreq > Config.lowerfreq)
          {
            const double lowfilfreq = Config.lowerfreq * Radial.GetDt();
            const double upfilfreq = Config.upperfreq * Radial.GetDt();
            SimpleBp Bandpass(lowfilfreq,upfilfreq);
            std::transform(Radial.GetData().begin(),Radial.GetData().end(),Radial.GetData().begin(),Bandpass);
            std::transform(Vertical.GetData().begin(),Vertical.GetData().end(),Vertical.GetData().begin(),Bandpass);
          }
        //initialize the object to calculate receiver functions with some essential parameters
        RecCalc
            Receiver(Config.shift, Config.sigma, Config.cc, false, rfmethod);
        //do we want the initial correlation peak to have an amplitude of 1 ?
        Receiver.SetNormalize(Config.normalize);
        SeismicDataComp RecFunc;
        //calculate the receiver function
        Receiver.CalcRecData(Radial, Vertical, RecFunc);
        //write out a sac file with the receiver function
        RecFunc.WriteAsSac(outfilename);
      } catch (FatalException &e)
      {
        cerr << e.what() << endl; // if something fails print error
        return -1; // and stop execution
      }
  }