Sdisp96Model.cpp

Go to the documentation of this file.

#include "Sdisp96Model.h"
#include <iostream>
#include <fstream>
#include <cassert>
#include <iterator>
#include <iomanip>
#include <boost/filesystem.hpp>
#include "FatalException.h"
using namespace std;

namespace gplib
  {
    Sdisp96Model::Sdisp96Model() :
      isSpherical(true)
      {
      }

    Sdisp96Model::~Sdisp96Model()
      {
      }

    Sdisp96Model& Sdisp96Model::operator=(const Sdisp96Model& source)
      {
        if (this == &source)
          return *this;
        SurfaceWaveModel::operator=(source);
        isSpherical = source.isSpherical;
        return *this;
      }

    Sdisp96Model::Sdisp96Model(const Sdisp96Model &Old) :
      SurfaceWaveModel(Old), isSpherical(Old.isSpherical)
      {

      }

    void Sdisp96Model::ReadModel(const std::string &filename)
      {

      }

    void Sdisp96Model::WriteModel(const std::string &filename) const
      {
        CheckConsistency();
        ofstream outfile(filename.c_str());
        outfile << "MODEL" << endl;
        outfile << GetName() << endl;
        outfile << "ISOTROPIC" << endl;
        outfile << "KGS" << endl;
        if (isSpherical)
          {
            outfile << "SPHERICAL EARTH" << endl;
          }
        else
          {
            outfile << "FLAT EARTH" << endl;
          }
        outfile << "1-D" << endl;
        outfile << "CONSTANT VELOCITY" << endl;
        outfile << "LINE08" << endl;
        outfile << "LINE09" << endl;
        outfile << "LINE10" << endl;
        outfile << "LINE11" << endl;
        outfile
            << "HR         VP         VS     RHO    QP    QS   ETAP    ETAS    FREFP    FREFS"
            << endl;
        const int nlayers = GetPvVelocities().size();
        for (int i = 0; i < nlayers; ++i)
          {
            outfile << setw(7) << setprecision(2) << GetThicknesses().at(i);
            outfile << setw(7) << setprecision(3) << GetPvVelocities().at(i);
            outfile << setw(7) << setprecision(3) << GetSvVelocities().at(i);
            outfile << setw(7) << setprecision(3) << GetDensities().at(i);
            outfile << setw(7) << setprecision(3) << GetQkappa().at(i);
            outfile << setw(7) << setprecision(3) << GetQmu().at(i);
            outfile << setw(7) << setprecision(3) << GetEta().at(i);
            outfile << setw(7) << setprecision(3) << GetEta().at(i);
            outfile << " 1.0  1.0" << endl;
          }
        if (!outfile.good())
          throw FatalException("Problem writing to file: " + filename);
      }

    void Sdisp96Model::WriteRunFile(const std::string &filename,
        const std::vector<double> periods) const
      {
        CheckConsistency();
        assert(periods.size() > 0);
        ofstream runfile(filename.c_str());
        ofstream periodfile((filename + ".per").c_str());
        copy(periods.begin(), periods.end(), ostream_iterator<double> (
            periodfile, "\n"));
        boost::filesystem::path Path(filename);
        string dirname(filename + "_dir");
        runfile << "#!/bin/bash" << endl;
        runfile << "mkdir " << dirname << endl;
        runfile << "mv " << filename << ".* " << dirname << endl;
        runfile << "cd " << dirname << endl;
        runfile << "sprep96  -M " << Path.leaf() << ".mod" << " -R -PARR "
            << Path.leaf() << ".per -NMOD 1 2>&1> /dev/null" << endl;
        runfile << "sdisp96 2>&1> /dev/null" << endl;
        runfile << "sdpsrf96 -R -ASC 2>&1> /dev/null" << endl;
        runfile << "cp SDISPR.ASC ../" << filename + ".asc96" << endl;

        if (!runfile.good())
          throw FatalException("Problem writing to file: " + filename);
      }
  }