MultiRecCalc.cpp

Go to the documentation of this file.

#include "MultiRecCalc.h"
#include "TimeFrequency.h"
#include "CFatalException.h"
#include "CRFSpecElement.h"
#include "TsSpectrum.h"
#include "WFunc.h"
#include "VecMat.h"
#include "Adaptors.h"

namespace ublas = boost::numeric::ublas;
using namespace std;
MultiRecCalc::MultiRecCalc(const int myshift, const double myomega,
    const double mysigma, const double myc) :
  c(myc), omega(myomega), sigma(mysigma), shift(myshift)
  {
  }

MultiRecCalc::~MultiRecCalc()
  {
  }


void MultiRecCalc::CalcRecData(const std::vector<SeismicDataComp> &RadComp,
    const std::vector<SeismicDataComp> &VerComp, SeismicDataComp &Receiver)
  {
    if (RadComp.size() != VerComp.size())
      throw CFatalException("Not the same number of stations for vertical and radial components !");
    const unsigned int compsize = RadComp.front().GetData().size();
    const unsigned int ncomp = RadComp.size();
    const double omegastep = 2. * PI/(RadComp.front().GetDt() * RadComp.front().GetData().size());
    gplib::rvec ConcatRad(compsize*ncomp), ConcatVer(compsize*ncomp);
    for (unsigned int i = 0; i < ncomp; ++i)
      {
        if (RadComp.at(i).GetData().size() != compsize || VerComp.at(i).GetData().size() != compsize)
          throw CFatalException("Incompatible sizes for Vertical or Radial component.");
        copy(RadComp.at(i).GetData().begin(), RadComp.at(i).GetData().end(), ConcatRad.begin()+i*compsize);
        copy(VerComp.at(i).GetData().begin(), VerComp.at(i).GetData().end(), ConcatVer.begin()+i*compsize);
      }
    gplib::cmat radtfspectrum = TimeFrequency(ConcatRad.begin(),
        ConcatRad.end(), compsize, Steep());
    gplib::cmat vertfspectrum = TimeFrequency(ConcatVer.begin(),
        ConcatVer.end(), compsize, Steep());
    gplib::cvec RecSpec(compsize/2+1);
    RecSpec *= 0.0;
    for (unsigned int i = 0; i < radtfspectrum.size1(); ++i)
      {
        complex<double> delta = *max_element(column(vertfspectrum,i).begin(), column(vertfspectrum,i).end(), gplib::absLess<tcomp,tcomp>()) *c;
        delta *= conj(delta);
        //delta /= radtfspectrum.size1();
        for (unsigned int j = 0; j < radtfspectrum.size2(); ++j)
          {
            RecSpec(j) += CalcSpectralElement(delta)(radtfspectrum(i, j),
                vertfspectrum(i, j));

          }
      }
    const double denom = 1./(4.0 * sigma*sigma);
    for (unsigned int i = 0; i < RecSpec.size(); ++i)
      {
        double freq = omegastep * i;
        RecSpec(i) *= exp(- (freq*freq)*denom);
      }
    Receiver.SetB(0.0);
    Receiver.SetDt(RadComp.front().GetDt());
    Receiver.GetData().assign(compsize, 0.0);
    TsSpectrum Spectrum(false);
    Spectrum.CalcTimeSeries(RecSpec.begin(), RecSpec.end(), Receiver.GetData().begin(), Receiver.GetData().end());
    const int shiftpoints =boost::numeric_cast<size_t>(shift/RadComp.front().GetDt());
    rotate(Receiver.GetData().begin(),Receiver.GetData().end()-shiftpoints,Receiver.GetData().end());
    const double maxamp = 1./ *max_element(Receiver.GetData().begin(),Receiver.GetData().end());
    transform(Receiver.GetData().begin(),Receiver.GetData().end(),Receiver.GetData().begin(),boost::bind(multiplies<double>(),_1,maxamp));
  }

---