AnisoSurfaceWaveModel.cpp

Go to the documentation of this file.

#include "AnisoSurfaceWaveModel.h"
#include <fstream>
#include <cassert>
#include <algorithm>
#include <numeric>
#include <boost/bind.hpp>
#include <string>
#include <boost/filesystem.hpp>
#include "FatalException.h"
#include <iostream>

namespace gplib
  {
    AnisoSurfaceWaveModel::AnisoSurfaceWaveModel()
      {
      }

    AnisoSurfaceWaveModel::~AnisoSurfaceWaveModel()
      {
      }

    void AnisoSurfaceWaveModel::ReadModel(const std::string &filename)
      {
        std::ifstream infile(filename.c_str());
        size_t nlayers = 0;
        const size_t titlelength = 512;
        char title[titlelength];
        //the first line of the model is the title that we don't care about
        infile.getline(title, titlelength);
        //the number of layers in the model does not include the bottom half-space
        infile >> nlayers;
        //so for us it is one layer more
        nlayers += 1;
        //there can be comments at the end of each line
        //so read the rest of the line, if nothing is there, it only reads the newline
        infile.getline(title, titlelength);
        //reserve memory for all the parameters
        thicknesses.assign(nlayers, 0.0);
        vp.assign(nlayers, 0.0);
        vs.assign(nlayers, 0.0);
        B.assign(nlayers, 0.0);
        C.assign(nlayers, 0.0);
        E.assign(nlayers, 0.0);
        theta.assign(nlayers, 0.0);
        phi.assign(nlayers, 0.0);
        densities.assign(nlayers, 0.0);
        for (size_t i = 0; i < nlayers; ++i)
          {
            //first line for each layer are the anisotropy angles
            infile >> theta.at(i) >> phi.at(i);
            //plus a comment
            infile.getline(title, titlelength);
            //then the other layer parameters
            infile >> thicknesses.at(i) >> vp.at(i) >> B.at(i) >> C.at(i)
                >> vs.at(i) >> E.at(i) >> densities.at(i);
            //transform from m to km for thickness and velocities and kg/m^3 to g/cm^3 for density
            thicknesses.at(i) /= 1000.0;
            vp.at(i) /= 1000.0;
            vs.at(i) /= 1000.0;
            densities.at(i) /= 1000.0;
            //swallow comment
            infile.getline(title, titlelength);
          }
        //the file contains depths we need thicknesses
        std::adjacent_difference(thicknesses.begin(), thicknesses.end(),
            thicknesses.begin());
      }

    void AnisoSurfaceWaveModel::WriteModel(const std::string &filename) const
      {
        const size_t nlayers = thicknesses.size();
        //check that all the layers have all parameters
        assert(nlayers> 0);
        assert(nlayers == vp.size());
        assert(nlayers == vs.size());
        assert(nlayers == B.size());
        assert(nlayers == C.size());
        assert(nlayers == E.size());
        assert(nlayers == theta.size());
        assert(nlayers == phi.size());
        assert(nlayers == densities.size());

        std::ofstream outfile(filename.c_str());
        //we have to write some title
        outfile << "Surface Wave model generated by gplib" << std::endl;
        //the code needs the number of layers without the half-space, we count the half-space
        outfile << nlayers - 1 << std::endl;
        //we have to transform thicknesses to depths
        double currdepth = 0.0;
        for (size_t i = 0; i < nlayers; ++i)
          {
            //write the anisotropy angles, the new line is important because the original file format can have comments
            outfile << theta.at(i) << " " << phi.at(i) << "\n";
            //calculate the depth to the bottom in m
            currdepth += thicknesses.at(i) * 1000.0;
            //write out the other parameters, again the newline matters
            //seismic anisotropic coefficient have to be positive (absolute value)
            outfile << currdepth << " " << vp.at(i) * 1000.0 << " " << fabs(
                B.at(i)) << " " << C.at(i) << " " << vs.at(i) * 1000.0 << " "
                << fabs(E.at(i)) << " " << densities.at(i) * 1000.0 << "\n";
          }
        if (!outfile.good())
          throw FatalException("Problem writing to file: " + filename);
      }

    void AnisoSurfaceWaveModel::WriteRunFile(const std::string &filename) const
      {
        std::ofstream runfile(filename.c_str());
        boost::filesystem::path Path(filename);
        std::string dirname(filename + "_dir");
        runfile << "#!/bin/bash" << std::endl;
        runfile << "mkdir " << dirname << std::endl;
        runfile << "mv " << filename << ".* " << dirname << std::endl;
        runfile << "cd " << dirname << std::endl;
        runfile << "cp " << Path.leaf() << ".mod animodel" << std::endl;
        runfile << "aniprop 2>&1> /dev/null" << std::endl;
        runfile << "cp out_cvel_R ../" << filename + ".cvel" << std::endl;
        if (!runfile.good())
          throw FatalException("Problem writing to file: " + filename);
      }

    void AnisoSurfaceWaveModel::WritePlot(const std::string &filename) const
      {
        std::ofstream outfile(filename.c_str());
        double cumthick = 0;
        const size_t nlayers = thicknesses.size();
        for (unsigned int i = 0; i < nlayers; ++i)
          {
            //we write S-wave velocity and one anisotropy parameter and strike
            //once for the top
            outfile << cumthick << "  " << vs.at(i) << "  " << B.at(i) << "  "
                << phi.at(i) << std::endl;
            cumthick += thicknesses.at(i);
            //and once for the bottom
            outfile << cumthick << " " << vs.at(i) << "  " << B.at(i) << "  "
                << phi.at(i) << std::endl;
          }
        //make sure we get a thick layer for the bottom half-space for plotting
        outfile << cumthick * 2.0 << " " << vs.back() << "  " << B.back()
            << "  " << phi.back() << std::endl;
      }
  }