RFVelCalc.cpp

Go to the documentation of this file.

#include "RFVelCalc.h"
#include "WFunc.h"
#include <numeric>
#include <fstream>
#include <cassert>
#include "CFatalException.h"
#include <gsl/gsl_math.h>

RFVelCalc::RFVelCalc(const double myomega, const double mysigma,
    const double myc, const RecCalc::trfmethod themethod) :
  RFCalculator(0, myomega, mysigma, myc, false, themethod)
  {
    RFCalculator.SetNormalize(false);
  }

RFVelCalc::~RFVelCalc()
  {
  }

RFVelCalc::RFVelCalc(const RFVelCalc &Old) :
  RFCalculator(Old.RFCalculator), Velocities(Old.Velocities),
      Periods(Old.Periods)
  {

  }

RFVelCalc& RFVelCalc::operator=(const RFVelCalc& source)
  {
    if (this == &source) return *this;
    RFCalculator = source.RFCalculator;
    Velocities = source.Velocities;
    Periods = source.Periods;
    return *this;
  }

void RFVelCalc::AbsVelCalc(const double slowness,
    const SeismicDataComp &RadRec, const SeismicDataComp &VerRec,
    ttsdata &AppVel)
  {
    ttsdata FilteredR(RadRec.GetData().size()), FilteredZ(VerRec.GetData().size());
    AppVel.clear();
    Velocities.clear();
    Periods.clear();
    assert(FilteredR.size() == FilteredZ.size());
    assert(gsl_fcmp(RadRec.GetDt(),VerRec.GetDt(),GSL_DBL_EPSILON)==0);
    for (double width = 0.1; width < 40; width *= 1.1) // we specify the width in seconds
      {
        double sigma = width/(RadRec.GetDt()*RadRec.GetData().size());
        ApplyWindow(RadRec.GetData().begin(), RadRec.GetData().end(), FilteredR.begin(), TruncCosSq(sigma));
        ApplyWindow(VerRec.GetData().begin(), VerRec.GetData().end(), FilteredZ.begin(), TruncCosSq(sigma));
        double rvalue = accumulate(FilteredR.begin(), FilteredR.end(), 0.0);
        double zvalue = accumulate(FilteredZ.begin(), FilteredZ.end(), 0.0);
        if (zvalue == 0.0)
          throw CFatalException("Zero amplitude for Z-Receiver function !");
        double angle = atan(rvalue/zvalue);
        double vsapp = sin(0.5 * angle)/slowness;
        AppVel.push_back(vsapp);
        Velocities.push_back(vsapp);
        Periods.push_back(width);
      }
  }

void RFVelCalc::CalcRFVelFromRec(const double slowness,
    const SeismicDataComp &RRec, const SeismicDataComp &VerComp,
    ttsdata &AppVel)
  {
    SeismicDataComp ZRecFunc;
    RFCalculator.CalcRecData(VerComp, VerComp, ZRecFunc);
    AbsVelCalc(slowness, RRec, ZRecFunc, AppVel);
  }

void RFVelCalc::CalcRFVel(const double slowness,
    const SeismicDataComp &RadComp, const SeismicDataComp &VerComp,
    ttsdata &AppVel)
  {
    SeismicDataComp RRecFunc, ZRecFunc;
    RFCalculator.CalcRecData(RadComp, VerComp, RRecFunc);
    RFCalculator.CalcRecData(VerComp, VerComp, ZRecFunc);
    AbsVelCalc(slowness, RRecFunc, ZRecFunc, AppVel);
  }

void RFVelCalc::WriteVelocities(const std::string filename)
  {
    std::ofstream outfile(filename.c_str());
    for (size_t i = 0; i < Velocities.size(); ++i)
      outfile << Periods.at(i) << " " << Velocities.at(i) << std::endl;
  }

---