calcrec.cpp

Go to the documentation of this file.

#include <iostream>
#include <string>
#include <algorithm>
#include <boost/math/constants/constants.hpp>
#include "Util.h"
#include "RecCalc.h"
#include "CalcRecConf.h"
#include "SeismicDataComp.h"
#include "FatalException.h"
#include "FilterFunc.h"
#include "statutils.h"
#include "miscfunc.h"
#include "NumUtil.h"
#include "RFVelCalc.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 1881 2014-02-25 16:20:58Z 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(argc, argv);
    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;

        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);
          }
        SeismicDataComp PComp(Vertical), SComp(Radial);
        if (Config.rotate)
          {
            RFVelCalc RFVel(Config.sigma, Config.cc, rfmethod);
            double p = Config.p;
            if (p < 0.0)
              {
                std::cerr << "Need to set ray parameter p " << std::endl;
                return 100;
              }
            ttsdata AppVel;
            //do the calculation and write to a file
            RFVel.CalcRFVel(p, Radial, Vertical, AppVel);
            double beta_est = Config.beta; //AppVel.at(0);
            double alpha_est = beta_est * sqrt(3.0);

            double qa = sqrt(std::pow(alpha_est, -2) - pow2(p));
            double qb = sqrt(std::pow(beta_est, -2) - pow2(p));
            double Vpz = (1.0 - 2.0 * pow2(beta_est) * pow2(p)) / (2.0 * alpha_est * qa);
            double Vsr = (1.0 - 2.0 * pow2(beta_est) * pow2(p)) / (2.0 * beta_est * qb);
            double Vpr = p * pow2(beta_est) / alpha_est;
            double Vsz = -p * beta_est;
            std::vector<double> VerNew(Vertical.GetData().size()), RadNew(
                Vertical.GetData().size());
            for (size_t i = 0; i < Vertical.GetData().size(); ++i)
              {
                VerNew.at(i) = Vpr * Radial.GetData().at(i)
                    + Vpz * Vertical.GetData().at(i);
                RadNew.at(i) = Vsr * Radial.GetData().at(i)
                    + Vsz * Vertical.GetData().at(i);
              }
            std::copy(RadNew.begin(), RadNew.end(), SComp.GetData().begin());
            std::copy(VerNew.begin(), VerNew.end(), PComp.GetData().begin());
            SComp.WriteAsSac("s.comp");
            PComp.WriteAsSac("p.comp");
          }
        //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
        if (Config.prec)
          {
            Receiver.CalcRecData(SComp, PComp, RecFunc);
          }
        else
          {
            Receiver.CalcRecData(PComp, SComp, 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
      }
  }