00001 #include "SurfaceWaveModel.h"
00002 #include <fstream>
00003 #include <cassert>
00004 #include <iterator>
00005 #include <iomanip>
00006 #include <numeric>
00007 #include <algorithm>
00008 #include <boost/bind.hpp>
00009 #include <boost/numeric/conversion/cast.hpp>
00010 #include <iostream>
00011 using namespace std;
00012
00013 namespace gplib
00014 {
00015 SurfaceWaveModel::SurfaceWaveModel() :
00016 refperiod(50), anisotropic(1), name("surf"), carddeckmodel(1),
00017 ninnercore(0), noutercore(0)
00018 {
00019 }
00020
00021 SurfaceWaveModel::~SurfaceWaveModel()
00022 {
00023 }
00024
00025 SurfaceWaveModel::SurfaceWaveModel(const SurfaceWaveModel &Old) :
00026 pvvelocities(Old.pvvelocities), phvelocities(Old.phvelocities),
00027 svvelocities(Old.svvelocities), shvelocities(Old.shvelocities),
00028 densities(Old.densities), thicknesses(Old.thicknesses), eta(Old.eta),
00029 qmu(Old.qmu), qkappa(Old.qkappa), refperiod(Old.refperiod),
00030 anisotropic(Old.anisotropic), name(Old.name), carddeckmodel(
00031 Old.carddeckmodel), ninnercore(Old.ninnercore), noutercore(
00032 Old.noutercore)
00033 {
00034
00035 }
00036
00037 SurfaceWaveModel& SurfaceWaveModel::operator=(
00038 const SurfaceWaveModel& source)
00039 {
00040 if (this == &source)
00041 return *this;
00042 source.CheckConsistency();
00043 pvvelocities.clear();
00044 phvelocities.clear();
00045 svvelocities.clear();
00046 shvelocities.clear();
00047 densities.clear();
00048 thicknesses.clear();
00049 eta.clear();
00050 qmu.clear();
00051 qkappa.clear();
00052 copy(source.pvvelocities.begin(), source.pvvelocities.end(),
00053 back_inserter(pvvelocities));
00054 copy(source.phvelocities.begin(), source.phvelocities.end(),
00055 back_inserter(phvelocities));
00056 copy(source.svvelocities.begin(), source.svvelocities.end(),
00057 back_inserter(svvelocities));
00058 copy(source.shvelocities.begin(), source.shvelocities.end(),
00059 back_inserter(shvelocities));
00060 copy(source.densities.begin(), source.densities.end(), back_inserter(
00061 densities));
00062 copy(source.thicknesses.begin(), source.thicknesses.end(),
00063 back_inserter(thicknesses));
00064 copy(source.eta.begin(), source.eta.end(), back_inserter(eta));
00065 copy(source.qmu.begin(), source.qmu.end(), back_inserter(qmu));
00066 copy(source.qkappa.begin(), source.qkappa.end(), back_inserter(qkappa));
00067
00068 refperiod = source.refperiod;
00069 anisotropic = source.anisotropic;
00070 name = source.name;
00071 carddeckmodel = source.carddeckmodel;
00072 ninnercore = source.ninnercore;
00073 noutercore = source.noutercore;
00074 return *this;
00075 }
00076
00077 void SurfaceWaveModel::CheckConsistency() const
00078 {
00079 assert(pvvelocities.size() == phvelocities.size());
00080 assert(pvvelocities.size() == svvelocities.size());
00081 assert(pvvelocities.size() == shvelocities.size());
00082 assert(pvvelocities.size() == densities.size());
00083 assert(pvvelocities.size() == thicknesses.size());
00084 assert(pvvelocities.size() == eta.size());
00085 assert(pvvelocities.size() == qmu.size());
00086 assert(pvvelocities.size() == qkappa.size());
00087 }
00088
00089 double SurfaceWaveModel::GetMaxDepth(const double depth)
00090 {
00091 assert(depth >= 0.0);
00092 if (depth < 10.0)
00093 return 1.0;
00094 if (depth < 20.0)
00095 return 2.0;
00096 if (depth < 100)
00097 return 5;
00098 if (depth < 400)
00099 return 10;
00100 return 100;
00101 }
00102
00103 void SurfaceWaveModel::AddDiscontinuity(const int index)
00104 {
00105 const unsigned int ninsert = 1;
00106 assert(index-1 < thicknesses.size());
00107 assert(index > 0);
00108
00109 thicknesses.insert(thicknesses.begin() + index, ninsert, 0.0);
00110 pvvelocities.insert(pvvelocities.begin() + index, ninsert,
00111 pvvelocities.at(index));
00112 phvelocities.insert(phvelocities.begin() + index, ninsert,
00113 phvelocities.at(index));
00114 svvelocities.insert(svvelocities.begin() + index, ninsert,
00115 svvelocities.at(index));
00116 shvelocities.insert(shvelocities.begin() + index, ninsert,
00117 shvelocities.at(index));
00118 densities.insert(densities.begin() + index, ninsert,
00119 densities.at(index));
00120 eta.insert(eta.begin() + index, ninsert, eta.at(index));
00121 qmu.insert(qmu.begin() + index, ninsert, qmu.at(index));
00122 qkappa.insert(qkappa.begin() + index, ninsert, qkappa.at(index));
00123 }
00124
00125 int SurfaceWaveModel::SplitLayer(const int index, const double maxthick)
00126 {
00127 assert(maxthick > 0.0);
00128 assert(index < thicknesses.size());
00129 double currthick = thicknesses.at(index);
00130 unsigned int ninsert = boost::numeric_cast<unsigned int>(floor(
00131 currthick / maxthick));
00132 if (ninsert > 0)
00133 {
00134 double diffthick = currthick - ninsert * maxthick;
00135 if (diffthick > 0.0)
00136
00137 {
00138 thicknesses.at(index) = diffthick;
00139 }
00140 else
00141
00142 {
00143 ninsert -= 1;
00144 thicknesses.at(index) = maxthick;
00145 }
00146 thicknesses.insert(thicknesses.begin() + index, ninsert, maxthick);
00147 pvvelocities.insert(pvvelocities.begin() + index, ninsert,
00148 pvvelocities.at(index));
00149 phvelocities.insert(phvelocities.begin() + index, ninsert,
00150 phvelocities.at(index));
00151 svvelocities.insert(svvelocities.begin() + index, ninsert,
00152 svvelocities.at(index));
00153 shvelocities.insert(shvelocities.begin() + index, ninsert,
00154 shvelocities.at(index));
00155 densities.insert(densities.begin() + index, ninsert, densities.at(
00156 index));
00157 eta.insert(eta.begin() + index, ninsert, eta.at(index));
00158 qmu.insert(qmu.begin() + index, ninsert, qmu.at(index));
00159 qkappa.insert(qkappa.begin() + index, ninsert, qkappa.at(index));
00160
00161 }
00162 if (index + ninsert + 1 < thicknesses.size())
00163 {
00164 AddDiscontinuity(index + ninsert + 1);
00165 return index + ninsert + 1;
00166 }
00167 return index + ninsert;
00168 }
00169
00170 void SurfaceWaveModel::MergeModel(const SurfaceWaveModel &Background)
00171 {
00172 CheckConsistency();
00173 Background.CheckConsistency();
00174 trealdata thisdepth(thicknesses.size(), 0.0), backgrounddepth(
00175 Background.thicknesses.size());
00176 partial_sum(thicknesses.begin(), thicknesses.end(), thisdepth.begin());
00177 partial_sum(Background.thicknesses.begin(),
00178 Background.thicknesses.end(), backgrounddepth.begin());
00179
00180 const unsigned int layerindex = distance(backgrounddepth.begin(),
00181 upper_bound(backgrounddepth.begin(), backgrounddepth.end(),
00182 thisdepth.back()));
00183
00184 assert(layerindex < Background.thicknesses.size());
00185 trealdata localthick(Background.thicknesses.begin() + layerindex,
00186 Background.thicknesses.end());
00187 localthick.front() = backgrounddepth.at(layerindex) - thisdepth.back();
00188
00189
00190 copy(localthick.begin(), localthick.end(), back_inserter(thicknesses));
00191 copy(Background.pvvelocities.begin() + layerindex,
00192 Background.pvvelocities.end(), back_inserter(pvvelocities));
00193 copy(Background.phvelocities.begin() + layerindex,
00194 Background.phvelocities.end(), back_inserter(phvelocities));
00195 copy(Background.svvelocities.begin() + layerindex,
00196 Background.svvelocities.end(), back_inserter(svvelocities));
00197 copy(Background.shvelocities.begin() + layerindex,
00198 Background.shvelocities.end(), back_inserter(shvelocities));
00199 copy(Background.densities.begin() + layerindex,
00200 Background.densities.end(), back_inserter(densities));
00201 copy(Background.eta.begin() + layerindex, Background.eta.end(),
00202 back_inserter(eta));
00203 copy(Background.qmu.begin() + layerindex, Background.qmu.end(),
00204 back_inserter(qmu));
00205 copy(Background.qkappa.begin() + layerindex, Background.qkappa.end(),
00206 back_inserter(qkappa));
00207 ninnercore = Background.ninnercore;
00208 noutercore = Background.noutercore;
00209 }
00210
00211 void SurfaceWaveModel::WriteLayer(std::ostream &stream,
00212 const unsigned int index, const trealdata &depth) const
00213 {
00214 stream << setw(10) << setprecision(6) << (6371 - depth.at(index))
00215 * 1000;
00216 stream << setw(10) << setprecision(6) << densities.at(index) * 1000
00217 << " " << pvvelocities.at(index) * 1000 << " ";
00218 stream << svvelocities.at(index) * 1000 << " " << qkappa.at(index)
00219 << " " << qmu.at(index) << " ";
00220 stream << phvelocities.at(index) * 1000 << " "
00221 << shvelocities.at(index) * 1000 << " " << eta.at(index) << " ";
00222 stream << depth.at(index);
00223 stream << endl;
00224 }
00225
00226 void SurfaceWaveModel::WritePlot(const std::string &filename) const
00227 {
00228 CheckConsistency();
00229 ofstream outfile(filename.c_str());
00230 double cumthick = 0;
00231 const unsigned int size = thicknesses.size();
00232 for (unsigned int i = 0; i < size; ++i)
00233 {
00234 outfile << cumthick << " " << svvelocities.at(i) << endl;
00235 cumthick += thicknesses.at(i);
00236 outfile << cumthick << " " << svvelocities.at(i) << endl;
00237 }
00238 outfile << cumthick * 2.0 << " " << svvelocities.at(size - 1) << endl;
00239 }
00240 }