SeismicModel.cpp

Go to the documentation of this file.

#include "SeismicModel.h"
#include <algorithm>
#include <fstream>
#include <iostream>
#include <numeric>
#include <cassert>

using namespace std;

namespace gplib
  {

    SeismicModel::SeismicModel(const int nlayers) :
      SourceDepth(0.0), dt(0.0), npts(0)
      {
        Init(nlayers);
      }

    SeismicModel::~SeismicModel()
      {
      }

    SeismicModel::SeismicModel(const SeismicModel& source)
      {
        //when we use the copy constructor the vectors
        //of the new object are empty so we can use a back_inserter
        copy(source.PVelocity.begin(), source.PVelocity.end(), back_inserter(
            this->PVelocity));
        copy(source.SVelocity.begin(), source.SVelocity.end(), back_inserter(
            this->SVelocity));
        copy(source.Density.begin(), source.Density.end(), back_inserter(
            this->Density));
        copy(source.Thickness.begin(), source.Thickness.end(), back_inserter(
            this->Thickness));
        copy(source.Qp.begin(), source.Qp.end(), back_inserter(this->Qp));
        copy(source.Qs.begin(), source.Qs.end(), back_inserter(this->Qs));
        SourceDepth = source.SourceDepth;
        dt = source.dt;
        npts = source.npts;
      }

    SeismicModel& SeismicModel::operator=(const SeismicModel& source)
      {
        if (this == &source)
          return *this;
        const int nelements = source.PVelocity.size();
        //for the regular copy operator we have to delete old data
        //so we allocate new space and copy the elements
        PVelocity.assign(nelements, 0);
        SVelocity.assign(nelements, 0);
        Density.assign(nelements, 0);
        Thickness.assign(nelements, 0);
        Qp.assign(nelements, 0);
        Qs.assign(nelements, 0);

        copy(source.PVelocity.begin(), source.PVelocity.end(),
            PVelocity.begin());
        copy(source.SVelocity.begin(), source.SVelocity.end(),
            SVelocity.begin());
        copy(source.Density.begin(), source.Density.end(), Density.begin());
        copy(source.Thickness.begin(), source.Thickness.end(),
            Thickness.begin());
        copy(source.Qp.begin(), source.Qp.end(), Qp.begin());
        copy(source.Qs.begin(), source.Qs.end(), Qs.begin());
        SourceDepth = source.SourceDepth;
        dt = source.dt;
        npts = source.npts;
        return *this;
      }

    void SeismicModel::Init(const int nlayers)
      {
        //with this function we can easily allocate
        //elements for the different vectors
        PVelocity.assign(nlayers, 0);
        SVelocity.assign(nlayers, 0);
        Density.assign(nlayers, 0);
        Thickness.assign(nlayers, 0);
        Qp.assign(nlayers, 0);
        Qs.assign(nlayers, 0);
      }
    //! Returns true if elem1 and elem2 agree within tolarance
    bool SeismicModel::FuzzComp(const double elem1, const double elem2,
        const double tolerance)
      {
        return (std::abs(elem1 - elem2) < tolerance);
      }

    int SeismicModel::FindLayer(const double depth)
      {
        //make sure the depth is positive
        assert(depth > 0.0);
        double currentbottom = 0;
        size_t i = 0;
        //we have to check that the bottom of the current layer
        //is deeper then what we are looking for
        //and that we don't look past the last index
        while (depth > currentbottom && i < Thickness.size())
          {
            currentbottom += Thickness.at(i);
            ++i;
          }
        //if we reach the end it means we are in the lower most
        //halfspace and return the index of the last layer
        return min(i, Thickness.size() - 1);
      }

    double SeismicModel::MatchSlowness(const double slowness,
        const tArrivalType mode)
      {
        double x = 0;
        double currdepth = 0;
        trealdata *velocity;
        const int sourcelayer = FindLayer(SourceDepth) - 1;

        if (mode == DirectS)
          velocity = &SVelocity;
        else
          velocity = &PVelocity;
        for (int i = 0; i < sourcelayer; ++i)
          {
            x += Thickness.at(i) * tan(asin(slowness * velocity->at(i)));
            currdepth += Thickness.at(i);
          }
        x += (SourceDepth - currdepth) * tan(asin(slowness * velocity->at(
            sourcelayer)));
        return x;
      }

    double SeismicModel::CalcTravelTime(const tArrivalType mode,
        const double sdepth, const double rdepth, const double p)
      {
        int sourceindex = FindLayer(sdepth);
        int recindex = FindLayer(rdepth);
        const double hdistance = MatchSlowness(p, mode);
        double currdepth = accumulate(Thickness.begin(), Thickness.begin()
            + recindex, 0.0);
        trealdata *velocity;
        if (mode == DirectS)
          velocity = &SVelocity;
        else
          velocity = &PVelocity;
        double t = p * hdistance;
        double eta = sqrt(1 / pow(velocity->at(recindex), 2) - p * p);
        t += (currdepth - rdepth) * eta;
        for (int i = recindex; i < sourceindex; ++i)
          {
            currdepth += Thickness.at(i);
            eta = sqrt(1 / pow(velocity->at(i), 2) - p * p);
            t += eta * Thickness.at(i);
          }
        eta = sqrt(1 / pow(velocity->at(sourceindex), 2) - p * p);
        t += (sdepth - currdepth) * eta;
        return t;
      }

    double SeismicModel::CalcArrival(const tArrivalType mode,
        const double recdist)
      {
        const double tolerance = 0.001;
        double angle = PI / 4.;
        double distance = 0;
        trealdata *velocity;
        double p = 0;
        double currdepth = 0;
        int preclevel = 3;

        if (mode == DirectS)
          velocity = &SVelocity;
        else
          velocity = &PVelocity;
        while (!FuzzComp(recdist, distance, tolerance))
          {
            currdepth = 0;
            p = sin(angle) / velocity->at(0);
            distance = MatchSlowness(p, mode);
            if (distance > recdist)
              angle -= PI * pow(0.5, preclevel);
            else
              angle += PI * pow(0.5, preclevel);
            preclevel++;
          }
        return CalcTravelTime(mode, SourceDepth, 0.0, p);
      }

    void SeismicModel::PlotVelWithErrors(const std::string &filename)
      {
        bool haveErrors = false;
        if (!SVelErrors.empty() && !ThickErrors.empty())
          haveErrors = true;
        ofstream outfile(filename.c_str());
        double cumthick = 0;
        const unsigned int size = Thickness.size();
        for (unsigned int i = 0; i < size; ++i)
          {
            outfile << cumthick << "  " << SVelocity.at(i);
            if (haveErrors)
              outfile << " " << ThickErrors.at(i) << " " << SVelErrors.at(i);
            outfile << endl;
            cumthick += Thickness.at(i);
            outfile << cumthick << " " << SVelocity.at(i);
            if (haveErrors)
              outfile << " " << ThickErrors.at(i) << " " << SVelErrors.at(i);
            outfile << endl;
          }
      }

    void SeismicModel::WritePlot(const std::string &filename)
      {
        ofstream outfile((filename + ".plot").c_str());
        double cumthick = 0;
        const unsigned int size = Thickness.size();
        for (unsigned int i = 0; i < size; ++i)
          {
            outfile << cumthick << "  " << SVelocity.at(i) << "  "
                << PVelocity.at(i) << "  " << Density.at(i) << endl;
            cumthick += Thickness.at(i);
            outfile << cumthick << " " << SVelocity.at(i) << "  "
                << PVelocity.at(i) << "  " << Density.at(i) << endl;
          }
        outfile << cumthick * 2.0 << " " << SVelocity.at(size - 1) << "  "
            << PVelocity.at(size - 1) << "  " << Density.at(size - 1) << endl;
      }
  }