SurfaceWaveModel.cpp

Go to the documentation of this file.
00001 #include "SurfaceWaveModel.h"
00002 #include <fstream>
00003 #include "CFatalException.h"
00004 #include <cassert>
00005 #include <iterator>
00006 #include <iomanip>
00007 #include <numeric>
00008 #include <algorithm>
00009 #include <boost/bind.hpp>
00010 #include <boost/numeric/conversion/cast.hpp>
00011 #include <iostream>
00012 using namespace std;
00013 
00014 SurfaceWaveModel::SurfaceWaveModel() :
00015   refperiod(50), anisotropic(1), name("surf"), carddeckmodel(1), ninnercore(0),
00016       noutercore(0)
00017   {
00018   }
00019 
00020 SurfaceWaveModel::~SurfaceWaveModel()
00021   {
00022   }
00023 
00024 SurfaceWaveModel::SurfaceWaveModel(const SurfaceWaveModel &Old) :
00025   pvvelocities(Old.pvvelocities), phvelocities(Old.phvelocities),
00026       svvelocities(Old.svvelocities), shvelocities(Old.shvelocities),
00027       densities(Old.densities), thicknesses(Old.thicknesses), eta(Old.eta),
00028       qmu(Old.qmu), qkappa(Old.qkappa), refperiod(Old.refperiod),
00029       anisotropic(Old.anisotropic), name(Old.name),
00030       carddeckmodel(Old.carddeckmodel), ninnercore(Old.ninnercore),
00031       noutercore(Old.noutercore)
00032   {
00033 
00034   }
00035 
00036 SurfaceWaveModel& SurfaceWaveModel::operator=(const SurfaceWaveModel& source)
00037   {
00038     if (this == &source) return *this;
00039     source.CheckConsistency();
00040     pvvelocities.clear();
00041     phvelocities.clear();
00042     svvelocities.clear();
00043     shvelocities.clear();
00044     densities.clear();
00045     thicknesses.clear();
00046     eta.clear();
00047     qmu.clear();
00048     qkappa.clear();
00049     copy(source.pvvelocities.begin(),source.pvvelocities.end(),back_inserter(pvvelocities));
00050     copy(source.phvelocities.begin(),source.phvelocities.end(),back_inserter(phvelocities));
00051     copy(source.svvelocities.begin(),source.svvelocities.end(),back_inserter(svvelocities));
00052     copy(source.shvelocities.begin(),source.shvelocities.end(),back_inserter(shvelocities));
00053     copy(source.densities.begin(),source.densities.end(),back_inserter(densities));
00054     copy(source.thicknesses.begin(),source.thicknesses.end(),back_inserter(thicknesses));
00055     copy(source.eta.begin(),source.eta.end(),back_inserter(eta));
00056     copy(source.qmu.begin(),source.qmu.end(),back_inserter(qmu));
00057     copy(source.qkappa.begin(),source.qkappa.end(),back_inserter(qkappa));
00058 
00059     refperiod = source.refperiod;
00060     anisotropic = source.anisotropic;
00061     name = source.name;
00062     carddeckmodel = source.carddeckmodel;
00063     ninnercore = source.ninnercore;
00064     noutercore = source.noutercore;
00065     return *this;
00066   }
00067 
00068 void SurfaceWaveModel::CheckConsistency() const
00069   {
00070     assert(pvvelocities.size() == phvelocities.size());
00071     assert(pvvelocities.size() == svvelocities.size());
00072     assert(pvvelocities.size() == shvelocities.size());
00073     assert(pvvelocities.size() == densities.size());
00074     assert(pvvelocities.size() == thicknesses.size());
00075     assert(pvvelocities.size() == eta.size());
00076     assert(pvvelocities.size() == qmu.size());
00077     assert(pvvelocities.size() == qkappa.size());
00078   }
00079 
00080 double SurfaceWaveModel::GetMaxDepth(const double depth)
00081   {
00082     assert(depth >= 0.0);
00083     if (depth < 10.0)
00084       return 1.0;
00085     if (depth < 20.0)
00086       return 2.0;
00087     if (depth < 100)
00088       return 5;
00089     if (depth < 400)
00090       return 10;
00091     return 100;
00092   }
00093 
00094 void SurfaceWaveModel::AddDiscontinuity(const int index)
00095   {
00096     const unsigned int ninsert = 1;
00097     assert(index-1 < thicknesses.size());
00098     assert(index > 0);
00099     
00100     thicknesses.insert(thicknesses.begin()+index, ninsert, 0.0);
00101     pvvelocities.insert(pvvelocities.begin()+index, ninsert,
00102         pvvelocities.at(index));
00103     phvelocities.insert(phvelocities.begin()+index, ninsert,
00104         phvelocities.at(index));
00105     svvelocities.insert(svvelocities.begin()+index, ninsert,
00106         svvelocities.at(index));
00107     shvelocities.insert(shvelocities.begin()+index, ninsert,
00108         shvelocities.at(index));
00109     densities.insert(densities.begin()+index, ninsert, densities.at(index));
00110     eta.insert(eta.begin()+index, ninsert, eta.at(index));
00111     qmu.insert(qmu.begin()+index, ninsert, qmu.at(index));
00112     qkappa.insert(qkappa.begin()+index, ninsert, qkappa.at(index));    
00113 }
00114 
00115 int SurfaceWaveModel::SplitLayer(const int index, const double maxthick)
00116   {
00117     assert(maxthick > 0.0);
00118     assert(index < thicknesses.size());
00119     double currthick = thicknesses.at(index);
00120     unsigned int ninsert = boost::numeric_cast<unsigned int>(floor(currthick/maxthick));
00121     if (ninsert > 0)
00122       {
00123         double diffthick = currthick - ninsert * maxthick;
00124         if (diffthick > 0.0) // we need diffthick extra layers and one with the difference
00125 
00126           {
00127             thicknesses.at(index) = diffthick;
00128           }
00129         else //they fit exactly
00130 
00131           {
00132             ninsert -=1;
00133             thicknesses.at(index) = maxthick;
00134           }
00135         thicknesses.insert(thicknesses.begin()+index,ninsert,maxthick);
00136         pvvelocities.insert(pvvelocities.begin()+index,ninsert,pvvelocities.at(index));
00137         phvelocities.insert(phvelocities.begin()+index,ninsert,phvelocities.at(index));
00138         svvelocities.insert(svvelocities.begin()+index,ninsert,svvelocities.at(index));
00139         shvelocities.insert(shvelocities.begin()+index,ninsert,shvelocities.at(index));
00140         densities.insert(densities.begin()+index,ninsert,densities.at(index));
00141         eta.insert(eta.begin()+index,ninsert,eta.at(index));
00142         qmu.insert(qmu.begin()+index,ninsert,qmu.at(index));
00143         qkappa.insert(qkappa.begin()+index,ninsert,qkappa.at(index));
00144 
00145       }
00146     if (index+ninsert+1 < thicknesses.size())
00147       {
00148         AddDiscontinuity(index+ninsert+1);
00149         return index+ninsert+1;
00150       }
00151     return index+ninsert;
00152   }
00153 
00154 void SurfaceWaveModel::MergeModel(const SurfaceWaveModel &Background)
00155   {
00156     CheckConsistency();
00157     Background.CheckConsistency();
00158     trealdata thisdepth(thicknesses.size(), 0.0),
00159         backgrounddepth(Background.thicknesses.size());
00160     partial_sum(thicknesses.begin(), thicknesses.end(), thisdepth.begin());
00161     partial_sum(Background.thicknesses.begin(), Background.thicknesses.end(),
00162         backgrounddepth.begin());
00163     //find the index of the layer in the backgroundmodel that lies below the last layer of this model  
00164     const unsigned int layerindex = distance(backgrounddepth.begin(),
00165         upper_bound(backgrounddepth.begin(), backgrounddepth.end(),
00166             thisdepth.back()));
00167     //cout << layerindex << " " << thisdepth.back() << " " << backgrounddepth.at(layerindex) << endl;
00168     assert(layerindex < Background.thicknesses.size());
00169     trealdata localthick(Background.thicknesses.begin()+layerindex,
00170         Background.thicknesses.end()); // we want to keep Background as it is
00171     localthick.front() = backgrounddepth.at(layerindex) - thisdepth.back(); //adjust thickness of bottom layer
00172     //copy(thicknesses.begin(),thicknesses.end(),ostream_iterator<double>(cout," "));
00173     //cout << endl;
00174     copy(localthick.begin(), localthick.end(), back_inserter(thicknesses));
00175     copy(Background.pvvelocities.begin()+layerindex,
00176         Background.pvvelocities.end(), back_inserter(pvvelocities));
00177     copy(Background.phvelocities.begin()+layerindex,
00178         Background.phvelocities.end(), back_inserter(phvelocities));
00179     copy(Background.svvelocities.begin()+layerindex,
00180         Background.svvelocities.end(), back_inserter(svvelocities));
00181     copy(Background.shvelocities.begin()+layerindex,
00182         Background.shvelocities.end(), back_inserter(shvelocities));
00183     copy(Background.densities.begin()+layerindex, Background.densities.end(),
00184         back_inserter(densities));
00185     copy(Background.eta.begin()+layerindex, Background.eta.end(),
00186         back_inserter(eta));
00187     copy(Background.qmu.begin()+layerindex, Background.qmu.end(),
00188         back_inserter(qmu));
00189     copy(Background.qkappa.begin()+layerindex, Background.qkappa.end(),
00190         back_inserter(qkappa));
00191     ninnercore = Background.ninnercore;
00192     noutercore = Background.noutercore;
00193   }
00194 
00195 void SurfaceWaveModel::ReadModel(const std::string &filename)
00196   {
00197     ifstream infile(filename.c_str());
00198     char line[1024];
00199     infile.getline(line, 1024);
00200     name = line;
00201     infile >> anisotropic >> refperiod >> carddeckmodel;
00202     if (!infile.good())
00203       throw CFatalException("Problem reading file: "+filename);
00204     unsigned int nlayers = 0;
00205     double dummy;
00206     infile >> nlayers >> ninnercore >> noutercore >> dummy >> dummy >> dummy
00207         >> dummy;
00208     unsigned int currlayer = 1;
00209     pvvelocities.assign(nlayers, 0.0);
00210     phvelocities.assign(nlayers, 0.0);
00211     svvelocities.assign(nlayers, 0.0);
00212     shvelocities.assign(nlayers, 0.0);
00213     densities.assign(nlayers, 0.0);
00214     thicknesses.assign(nlayers, 0.0);
00215     eta.assign(nlayers, 0.0);
00216     qmu.assign(nlayers, 0.0);
00217     qkappa.assign(nlayers, 0.0);
00218     trealdata radius(nlayers, 0.0);
00219     int index = nlayers-1;
00220     double currrad;
00221     while (infile.good() && index >= 0)
00222       {
00223         infile >> currrad;
00224         radius.at(index) = 6371 - currrad / 1000.0;
00225         if (infile.good())
00226           {
00227             infile >> densities.at(index) >> pvvelocities.at(index);
00228             infile >> svvelocities.at(index) >> qkappa.at(index)
00229                 >> qmu.at(index);
00230             infile >> phvelocities.at(index) >> shvelocities.at(index)
00231                 >> eta.at(index);
00232             infile >> dummy;
00233             currlayer++;
00234             index--;
00235           }
00236 
00237       }
00238     transform(densities.begin(), densities.end(), densities.begin(),
00239         boost::bind(multiplies<double>(), _1, 1./1000.0));
00240     transform(svvelocities.begin(), svvelocities.end(), svvelocities.begin(),
00241         boost::bind(multiplies<double>(), _1, 1./1000.0));
00242     transform(shvelocities.begin(), shvelocities.end(), shvelocities.begin(),
00243         boost::bind(multiplies<double>(), _1, 1./1000.0));
00244     transform(pvvelocities.begin(), pvvelocities.end(), pvvelocities.begin(),
00245         boost::bind(multiplies<double>(), _1, 1./1000.0));
00246     transform(phvelocities.begin(), phvelocities.end(), phvelocities.begin(),
00247         boost::bind(multiplies<double>(), _1, 1./1000.0));
00248 
00249     //if (currlayer != nlayers)
00250     //  throw CFatalException("Problem reading file: "+filename);
00251     adjacent_difference(radius.begin(), radius.end(), thicknesses.begin());
00252   }
00253 
00254 void SurfaceWaveModel::WriteLayer(std::ostream &stream,
00255     const unsigned int index, const trealdata &depth) const
00256   {
00257     stream << setw(10) << setprecision(6) << (6371 - depth.at(index)) * 1000;
00258     stream << setw(10) << setprecision(6) << densities.at(index) * 1000 << " "
00259         <<pvvelocities.at(index) * 1000 << " ";
00260     stream << svvelocities.at(index) * 1000 <<" " << qkappa.at(index) <<" "
00261         << qmu.at(index) << " ";
00262     stream << phvelocities.at(index) * 1000 <<" " << shvelocities.at(index)
00263         * 1000 <<" " << eta.at(index) << " ";
00264     stream << depth.at(index); // have to write depth here
00265     stream << endl;
00266   }
00267 
00268 void SurfaceWaveModel::WriteModel(const std::string &filename) const
00269   {
00270     CheckConsistency();
00271     const int nlayers = pvvelocities.size();
00272     ofstream outfile(filename.c_str());
00273 
00274     outfile << "anisotropic " << name << endl;
00275     outfile << anisotropic << " " << refperiod << " " << carddeckmodel << endl;
00276 
00277     trealdata depth(nlayers, 0.0);
00278     partial_sum(thicknesses.begin(), thicknesses.end(), depth.begin());
00279 
00280     outfile << nlayers << " " << ninnercore << " " << noutercore << " "
00281         << nlayers << " " << nlayers << " " << nlayers << " 0" << endl;
00282     for (int i = nlayers-1; i >= 1; --i)
00283       {
00284         WriteLayer(outfile, i, depth);
00285       }
00286     outfile << setw(10) << setprecision(6) << 6371 * 1000;
00287     outfile << setw(10) << setprecision(6) << densities.front() * 1000 << " "
00288         <<pvvelocities.front() * 1000 << " ";
00289     outfile << svvelocities.front() * 1000 <<" " << qkappa.front() <<" "
00290         << qmu.front() << " ";
00291     outfile << phvelocities.front() * 1000 <<" " << shvelocities.front() * 1000
00292         <<" " << eta.front() << " ";
00293     outfile << 0.0 << endl;
00294     if (!outfile.good())
00295       throw CFatalException("Problem writing to file: "+filename);
00296   }
00297 
00298 void SurfaceWaveModel::WritePlot(const std::string &filename) const
00299   {
00300     CheckConsistency();
00301     ofstream outfile(filename.c_str());
00302     double cumthick = 0;
00303     const unsigned int size = thicknesses.size();
00304     for (unsigned int i = 0; i < size; ++i)
00305       {
00306         outfile << cumthick << "  " << svvelocities.at(i) << endl;
00307         cumthick += thicknesses.at(i);
00308         outfile << cumthick << " " << svvelocities.at(i) << endl;
00309       }
00310     outfile << cumthick*2.0 << " " << svvelocities.at(size-1) << endl;
00311   }
00312 
00313 void SurfaceWaveModel::WriteRunFile(const std::string &filename,
00314     const std::vector<double> periods) const
00315   {
00316     assert(periods.size() > 0);
00317     ofstream runfile(filename.c_str());
00318     runfile << "#!/bin/bash" << endl;
00319     runfile << "rayleigh << eof > /dev/null" << endl;
00320     runfile << (filename+".mod").c_str() << endl;
00321     runfile << "out" << endl;
00322     runfile << (filename+".pvel").c_str() << endl;
00323     runfile << periods.size() << endl;
00324     copy(periods.begin(), periods.end(), ostream_iterator<double>(runfile, " "));
00325     runfile << endl << "eof" << endl;
00326     if (!runfile.good())
00327       throw CFatalException("Problem writing to file: "+filename);
00328   }
00329 

Generated on Fri Jul 4 15:30:21 2008 for GPLIB++ by  doxygen 1.5.5