00001 #include "Sdisp96Model.h"
00002 #include <iostream>
00003 #include <fstream>
00004 #include <cassert>
00005 #include <iterator>
00006 #include <iomanip>
00007 #include <boost/filesystem.hpp>
00008 #include "CFatalException.h"
00009 using namespace std;
00010
00011 Sdisp96Model::Sdisp96Model() :
00012 isSpherical(true)
00013 {
00014 }
00015
00016 Sdisp96Model::~Sdisp96Model()
00017 {
00018 }
00019
00020 Sdisp96Model& Sdisp96Model::operator=(const Sdisp96Model& source)
00021 {
00022 if (this == &source)
00023 return *this;
00024 SurfaceWaveModel::operator=(source);
00025 isSpherical = source.isSpherical;
00026 return *this;
00027 }
00028
00029 Sdisp96Model::Sdisp96Model(const Sdisp96Model &Old) :
00030 SurfaceWaveModel(Old), isSpherical(Old.isSpherical)
00031 {
00032
00033 }
00034
00035 void Sdisp96Model::ReadModel(const std::string &filename)
00036 {
00037
00038 }
00039
00040 void Sdisp96Model::WriteModel(const std::string &filename) const
00041 {
00042 CheckConsistency();
00043 ofstream outfile(filename.c_str());
00044 outfile << "MODEL" << endl;
00045 outfile << GetName() << endl;
00046 outfile << "ISOTROPIC" << endl;
00047 outfile << "KGS" << endl;
00048 if (isSpherical)
00049 {
00050 outfile << "SPHERICAL EARTH" << endl;
00051 }
00052 else
00053 {
00054 outfile << "FLAT EARTH" << endl;
00055 }
00056 outfile << "1-D" << endl;
00057 outfile << "CONSTANT VELOCITY" << endl;
00058 outfile << "LINE08" << endl;
00059 outfile << "LINE09" << endl;
00060 outfile << "LINE10" << endl;
00061 outfile << "LINE11" << endl;
00062 outfile
00063 << "HR VP VS RHO QP QS ETAP ETAS FREFP FREFS"
00064 << endl;
00065 const int nlayers = GetPvVelocities().size();
00066 for (int i = 0; i < nlayers; ++i)
00067 {
00068 outfile << setw(7) << setprecision(3) << GetThicknesses().at(i);
00069 outfile << setw(7) << setprecision(3) << GetPvVelocities().at(i);
00070 outfile << setw(7) << setprecision(3) << GetSvVelocities().at(i);
00071 outfile << setw(7) << setprecision(3) << GetDensities().at(i);
00072 outfile << setw(7) << setprecision(3) << GetQkappa().at(i);
00073 outfile << setw(7) << setprecision(3) << GetQmu().at(i);
00074 outfile << setw(7) << setprecision(3) << GetEta().at(i);
00075 outfile << setw(7) << setprecision(3) << GetEta().at(i);
00076 outfile << " 1.0 1.0" << endl;
00077 }
00078 }
00079
00080 void Sdisp96Model::WriteRunFile(const std::string &filename,
00081 const std::vector<double> periods) const
00082 {
00083 CheckConsistency();
00084 assert(periods.size() > 0);
00085 ofstream runfile(filename.c_str());
00086 ofstream periodfile((filename+".per").c_str());
00087 copy(periods.begin(), periods.end(), ostream_iterator<double>(periodfile,
00088 "\n"));
00089 boost::filesystem::path Path(filename);
00090 string dirname(filename+"_dir");
00091 runfile << "#!/bin/bash" << endl;
00092 runfile << "mkdir " << dirname << endl;
00093 runfile << "mv " << filename << ".* " << dirname << endl;
00094 runfile << "cd " << dirname << endl;
00095 runfile << "sprep96 -M "<< filename << "_dir/" << Path.leaf() << ".mod" << " -R -PARR "<<filename << "_dir/" << Path.leaf()
00096 <<".per -NMOD 1 2>&1> /dev/null" << endl;
00097 runfile << "sdisp96 2>&1> /dev/null" << endl;
00098 runfile << "sdpsrf96 -R -ASC 2>&1> /dev/null" << endl;
00099 runfile << "cp SDISPR.ASC " <<filename+".asc96" << endl;
00100
00101 if (!runfile.good())
00102 throw CFatalException("Problem writing to file: "+filename);
00103 }
