SurfaceWaveModel.cpp

Go to the documentation of this file.

#include "SurfaceWaveModel.h"
#include <fstream>
#include "CFatalException.h"
#include <cassert>
#include <iterator>
#include <iomanip>
#include <numeric>
#include <algorithm>
#include <boost/bind.hpp>
#include <boost/numeric/conversion/cast.hpp>
#include <iostream>
using namespace std;

SurfaceWaveModel::SurfaceWaveModel() :
  refperiod(50), anisotropic(1), name("surf"), carddeckmodel(1), ninnercore(0),
      noutercore(0)
  {
  }

SurfaceWaveModel::~SurfaceWaveModel()
  {
  }

SurfaceWaveModel::SurfaceWaveModel(const SurfaceWaveModel &Old) :
  pvvelocities(Old.pvvelocities), phvelocities(Old.phvelocities),
      svvelocities(Old.svvelocities), shvelocities(Old.shvelocities),
      densities(Old.densities), thicknesses(Old.thicknesses), eta(Old.eta),
      qmu(Old.qmu), qkappa(Old.qkappa), refperiod(Old.refperiod),
      anisotropic(Old.anisotropic), name(Old.name),
      carddeckmodel(Old.carddeckmodel), ninnercore(Old.ninnercore),
      noutercore(Old.noutercore)
  {

  }

SurfaceWaveModel& SurfaceWaveModel::operator=(const SurfaceWaveModel& source)
  {
    if (this == &source) return *this;
    source.CheckConsistency();
    pvvelocities.clear();
    phvelocities.clear();
    svvelocities.clear();
    shvelocities.clear();
    densities.clear();
    thicknesses.clear();
    eta.clear();
    qmu.clear();
    qkappa.clear();
    copy(source.pvvelocities.begin(),source.pvvelocities.end(),back_inserter(pvvelocities));
    copy(source.phvelocities.begin(),source.phvelocities.end(),back_inserter(phvelocities));
    copy(source.svvelocities.begin(),source.svvelocities.end(),back_inserter(svvelocities));
    copy(source.shvelocities.begin(),source.shvelocities.end(),back_inserter(shvelocities));
    copy(source.densities.begin(),source.densities.end(),back_inserter(densities));
    copy(source.thicknesses.begin(),source.thicknesses.end(),back_inserter(thicknesses));
    copy(source.eta.begin(),source.eta.end(),back_inserter(eta));
    copy(source.qmu.begin(),source.qmu.end(),back_inserter(qmu));
    copy(source.qkappa.begin(),source.qkappa.end(),back_inserter(qkappa));

    refperiod = source.refperiod;
    anisotropic = source.anisotropic;
    name = source.name;
    carddeckmodel = source.carddeckmodel;
    ninnercore = source.ninnercore;
    noutercore = source.noutercore;
    return *this;
  }

void SurfaceWaveModel::CheckConsistency() const
  {
    assert(pvvelocities.size() == phvelocities.size());
    assert(pvvelocities.size() == svvelocities.size());
    assert(pvvelocities.size() == shvelocities.size());
    assert(pvvelocities.size() == densities.size());
    assert(pvvelocities.size() == thicknesses.size());
    assert(pvvelocities.size() == eta.size());
    assert(pvvelocities.size() == qmu.size());
    assert(pvvelocities.size() == qkappa.size());
  }

double SurfaceWaveModel::GetMaxDepth(const double depth)
  {
    assert(depth >= 0.0);
    if (depth < 10.0)
      return 1.0;
    if (depth < 20.0)
      return 2.0;
    if (depth < 100)
      return 5;
    if (depth < 400)
      return 10;
    return 100;
  }

void SurfaceWaveModel::AddDiscontinuity(const int index)
  {
    const unsigned int ninsert = 1;
    assert(index-1 < thicknesses.size());
    assert(index > 0);
    
    thicknesses.insert(thicknesses.begin()+index, ninsert, 0.0);
    pvvelocities.insert(pvvelocities.begin()+index, ninsert,
        pvvelocities.at(index));
    phvelocities.insert(phvelocities.begin()+index, ninsert,
        phvelocities.at(index));
    svvelocities.insert(svvelocities.begin()+index, ninsert,
        svvelocities.at(index));
    shvelocities.insert(shvelocities.begin()+index, ninsert,
        shvelocities.at(index));
    densities.insert(densities.begin()+index, ninsert, densities.at(index));
    eta.insert(eta.begin()+index, ninsert, eta.at(index));
    qmu.insert(qmu.begin()+index, ninsert, qmu.at(index));
    qkappa.insert(qkappa.begin()+index, ninsert, qkappa.at(index));    
}

int SurfaceWaveModel::SplitLayer(const int index, const double maxthick)
  {
    assert(maxthick > 0.0);
    assert(index < thicknesses.size());
    double currthick = thicknesses.at(index);
    unsigned int ninsert = boost::numeric_cast<unsigned int>(floor(currthick/maxthick));
    if (ninsert > 0)
      {
        double diffthick = currthick - ninsert * maxthick;
        if (diffthick > 0.0) // we need diffthick extra layers and one with the difference

          {
            thicknesses.at(index) = diffthick;
          }
        else //they fit exactly

          {
            ninsert -=1;
            thicknesses.at(index) = maxthick;
          }
        thicknesses.insert(thicknesses.begin()+index,ninsert,maxthick);
        pvvelocities.insert(pvvelocities.begin()+index,ninsert,pvvelocities.at(index));
        phvelocities.insert(phvelocities.begin()+index,ninsert,phvelocities.at(index));
        svvelocities.insert(svvelocities.begin()+index,ninsert,svvelocities.at(index));
        shvelocities.insert(shvelocities.begin()+index,ninsert,shvelocities.at(index));
        densities.insert(densities.begin()+index,ninsert,densities.at(index));
        eta.insert(eta.begin()+index,ninsert,eta.at(index));
        qmu.insert(qmu.begin()+index,ninsert,qmu.at(index));
        qkappa.insert(qkappa.begin()+index,ninsert,qkappa.at(index));

      }
    if (index+ninsert+1 < thicknesses.size())
      {
        AddDiscontinuity(index+ninsert+1);
        return index+ninsert+1;
      }
    return index+ninsert;
  }

void SurfaceWaveModel::MergeModel(const SurfaceWaveModel &Background)
  {
    CheckConsistency();
    Background.CheckConsistency();
    trealdata thisdepth(thicknesses.size(), 0.0),
        backgrounddepth(Background.thicknesses.size());
    partial_sum(thicknesses.begin(), thicknesses.end(), thisdepth.begin());
    partial_sum(Background.thicknesses.begin(), Background.thicknesses.end(),
        backgrounddepth.begin());
    //find the index of the layer in the backgroundmodel that lies below the last layer of this model  
    const unsigned int layerindex = distance(backgrounddepth.begin(),
        upper_bound(backgrounddepth.begin(), backgrounddepth.end(),
            thisdepth.back()));
    //cout << layerindex << " " << thisdepth.back() << " " << backgrounddepth.at(layerindex) << endl;
    assert(layerindex < Background.thicknesses.size());
    trealdata localthick(Background.thicknesses.begin()+layerindex,
        Background.thicknesses.end()); // we want to keep Background as it is
    localthick.front() = backgrounddepth.at(layerindex) - thisdepth.back(); //adjust thickness of bottom layer
    //copy(thicknesses.begin(),thicknesses.end(),ostream_iterator<double>(cout," "));
    //cout << endl;
    copy(localthick.begin(), localthick.end(), back_inserter(thicknesses));
    copy(Background.pvvelocities.begin()+layerindex,
        Background.pvvelocities.end(), back_inserter(pvvelocities));
    copy(Background.phvelocities.begin()+layerindex,
        Background.phvelocities.end(), back_inserter(phvelocities));
    copy(Background.svvelocities.begin()+layerindex,
        Background.svvelocities.end(), back_inserter(svvelocities));
    copy(Background.shvelocities.begin()+layerindex,
        Background.shvelocities.end(), back_inserter(shvelocities));
    copy(Background.densities.begin()+layerindex, Background.densities.end(),
        back_inserter(densities));
    copy(Background.eta.begin()+layerindex, Background.eta.end(),
        back_inserter(eta));
    copy(Background.qmu.begin()+layerindex, Background.qmu.end(),
        back_inserter(qmu));
    copy(Background.qkappa.begin()+layerindex, Background.qkappa.end(),
        back_inserter(qkappa));
    ninnercore = Background.ninnercore;
    noutercore = Background.noutercore;
  }

void SurfaceWaveModel::ReadModel(const std::string &filename)
  {
    ifstream infile(filename.c_str());
    char line[1024];
    infile.getline(line, 1024);
    name = line;
    infile >> anisotropic >> refperiod >> carddeckmodel;
    if (!infile.good())
      throw CFatalException("Problem reading file: "+filename);
    unsigned int nlayers = 0;
    double dummy;
    infile >> nlayers >> ninnercore >> noutercore >> dummy >> dummy >> dummy
        >> dummy;
    unsigned int currlayer = 1;
    pvvelocities.assign(nlayers, 0.0);
    phvelocities.assign(nlayers, 0.0);
    svvelocities.assign(nlayers, 0.0);
    shvelocities.assign(nlayers, 0.0);
    densities.assign(nlayers, 0.0);
    thicknesses.assign(nlayers, 0.0);
    eta.assign(nlayers, 0.0);
    qmu.assign(nlayers, 0.0);
    qkappa.assign(nlayers, 0.0);
    trealdata radius(nlayers, 0.0);
    int index = nlayers-1;
    double currrad;
    while (infile.good() && index >= 0)
      {
        infile >> currrad;
        radius.at(index) = 6371 - currrad / 1000.0;
        if (infile.good())
          {
            infile >> densities.at(index) >> pvvelocities.at(index);
            infile >> svvelocities.at(index) >> qkappa.at(index)
                >> qmu.at(index);
            infile >> phvelocities.at(index) >> shvelocities.at(index)
                >> eta.at(index);
            infile >> dummy;
            currlayer++;
            index--;
          }

      }
    transform(densities.begin(), densities.end(), densities.begin(),
        boost::bind(multiplies<double>(), _1, 1./1000.0));
    transform(svvelocities.begin(), svvelocities.end(), svvelocities.begin(),
        boost::bind(multiplies<double>(), _1, 1./1000.0));
    transform(shvelocities.begin(), shvelocities.end(), shvelocities.begin(),
        boost::bind(multiplies<double>(), _1, 1./1000.0));
    transform(pvvelocities.begin(), pvvelocities.end(), pvvelocities.begin(),
        boost::bind(multiplies<double>(), _1, 1./1000.0));
    transform(phvelocities.begin(), phvelocities.end(), phvelocities.begin(),
        boost::bind(multiplies<double>(), _1, 1./1000.0));

    //if (currlayer != nlayers)
    //  throw CFatalException("Problem reading file: "+filename);
    adjacent_difference(radius.begin(), radius.end(), thicknesses.begin());
  }

void SurfaceWaveModel::WriteLayer(std::ostream &stream,
    const unsigned int index, const trealdata &depth) const
  {
    stream << setw(10) << setprecision(6) << (6371 - depth.at(index)) * 1000;
    stream << setw(10) << setprecision(6) << densities.at(index) * 1000 << " "
        <<pvvelocities.at(index) * 1000 << " ";
    stream << svvelocities.at(index) * 1000 <<" " << qkappa.at(index) <<" "
        << qmu.at(index) << " ";
    stream << phvelocities.at(index) * 1000 <<" " << shvelocities.at(index)
        * 1000 <<" " << eta.at(index) << " ";
    stream << depth.at(index); // have to write depth here
    stream << endl;
  }

void SurfaceWaveModel::WriteModel(const std::string &filename) const
  {
    CheckConsistency();
    const int nlayers = pvvelocities.size();
    ofstream outfile(filename.c_str());

    outfile << "anisotropic " << name << endl;
    outfile << anisotropic << " " << refperiod << " " << carddeckmodel << endl;

    trealdata depth(nlayers, 0.0);
    partial_sum(thicknesses.begin(), thicknesses.end(), depth.begin());

    outfile << nlayers << " " << ninnercore << " " << noutercore << " "
        << nlayers << " " << nlayers << " " << nlayers << " 0" << endl;
    for (int i = nlayers-1; i >= 1; --i)
      {
        WriteLayer(outfile, i, depth);
      }
    outfile << setw(10) << setprecision(6) << 6371 * 1000;
    outfile << setw(10) << setprecision(6) << densities.front() * 1000 << " "
        <<pvvelocities.front() * 1000 << " ";
    outfile << svvelocities.front() * 1000 <<" " << qkappa.front() <<" "
        << qmu.front() << " ";
    outfile << phvelocities.front() * 1000 <<" " << shvelocities.front() * 1000
        <<" " << eta.front() << " ";
    outfile << 0.0 << endl;
    if (!outfile.good())
      throw CFatalException("Problem writing to file: "+filename);
  }

void SurfaceWaveModel::WritePlot(const std::string &filename) const
  {
    CheckConsistency();
    ofstream outfile(filename.c_str());
    double cumthick = 0;
    const unsigned int size = thicknesses.size();
    for (unsigned int i = 0; i < size; ++i)
      {
        outfile << cumthick << "  " << svvelocities.at(i) << endl;
        cumthick += thicknesses.at(i);
        outfile << cumthick << " " << svvelocities.at(i) << endl;
      }
    outfile << cumthick*2.0 << " " << svvelocities.at(size-1) << endl;
  }

void SurfaceWaveModel::WriteRunFile(const std::string &filename,
    const std::vector<double> periods) const
  {
    assert(periods.size() > 0);
    ofstream runfile(filename.c_str());
    runfile << "#!/bin/bash" << endl;
    runfile << "rayleigh << eof > /dev/null" << endl;
    runfile << (filename+".mod").c_str() << endl;
    runfile << "out" << endl;
    runfile << (filename+".pvel").c_str() << endl;
    runfile << periods.size() << endl;
    copy(periods.begin(), periods.end(), ostream_iterator<double>(runfile, " "));
    runfile << endl << "eof" << endl;
    if (!runfile.good())
      throw CFatalException("Problem writing to file: "+filename);
  }

---