00001 #include "AnisoSurfaceWaveModel.h" 00002 #include <fstream> 00003 #include <cassert> 00004 #include <algorithm> 00005 #include <numeric> 00006 #include <boost/bind.hpp> 00007 #include <string> 00008 #include <boost/filesystem.hpp> 00009 #include "FatalException.h" 00010 #include <iostream> 00011 00012 namespace gplib 00013 { 00014 00015 AnisoSurfaceWaveModel::AnisoSurfaceWaveModel() 00016 { 00017 } 00018 00019 AnisoSurfaceWaveModel::~AnisoSurfaceWaveModel() 00020 { 00021 } 00022 00023 // Function Read Model if "gradient" in half-space 00024 /*void AnisoSurfaceWaveModel::ReadModel(const std::string &filename) 00025 { 00026 std::ifstream infile(filename.c_str()); 00027 size_t nlayers = 0; 00028 const size_t titlelength = 512; 00029 char title[titlelength]; 00030 //the first line of the model is the title that we don't care about 00031 infile.getline(title, titlelength); 00032 //the number of layers in the model does not include the bottom half-space 00033 infile >> nlayers; 00034 //so for us it is one layer more 00035 nlayers += 1; 00036 //there can be comments at the end of each line 00037 //so read the rest of the line, if nothing is there, it only reads the newline 00038 infile.getline(title, titlelength); 00039 //reserve memory for all the parameters 00040 00041 // remove comment for a "gradient" in half-space 00042 thicknesses.assign(nlayers -4, 0.0); 00043 vp.assign(nlayers -4, 0.0); 00044 vs.assign(nlayers -4, 0.0); 00045 B.assign(nlayers -4, 0.0); 00046 C.assign(nlayers -4, 0.0); 00047 E.assign(nlayers -4, 0.0); 00048 theta.assign(nlayers -4, 0.0); 00049 phi.assign(nlayers -4, 0.0); 00050 densities.assign(nlayers -4, 0.0); 00051 for (size_t i = 0; i < nlayers -5; ++i) 00052 { 00053 //first line for each layer are the anisotropy angles 00054 infile >> theta.at(i) >> phi.at(i); 00055 //plus a comment 00056 infile.getline(title, titlelength); 00057 //then the other layer parameters 00058 infile >> thicknesses.at(i) >> vp.at(i) >> B.at(i) >> C.at(i) >> vs.at( 00059 i) >> E.at(i) >> densities.at(i); 00060 //transform from m to km for thickness and velocities and kg/m^3 to g/cm^3 for density 00061 thicknesses.at(i) /= 1000.0; 00062 vp.at(i) /= 1000.0; 00063 vs.at(i) /= 1000.0; 00064 densities.at(i) /= 1000.0; 00065 //swallow comment 00066 infile.getline(title, titlelength); 00067 } 00068 // then we read the last layer divided in 5 layers of same thickness 00069 //for the layer n-5 00070 infile >> theta.at(nlayers - 5) >> phi.at(nlayers - 5); 00071 infile.getline(title, titlelength); 00072 infile >> thicknesses.at(nlayers - 5) >> vp.at(nlayers - 5) >> B.at(nlayers - 5) >> C.at(nlayers - 5) >> vs.at( 00073 nlayers - 5) >> E.at(nlayers - 5) >> densities.at(nlayers - 5); 00074 thicknesses.at(nlayers - 5) /= 1000.0; 00075 vp.at(nlayers - 5) /= 1000.0; 00076 vs.at(nlayers - 5) /= 1000.0; 00077 densities.at(nlayers - 5) /= 1000.0; 00078 infile.getline(title, titlelength); 00079 //for the layer n-4 00080 infile >> theta.at(nlayers - 5) >> phi.at(nlayers - 5); 00081 infile.getline(title, titlelength); 00082 infile >> thicknesses.at(nlayers - 5) >> vp.at(nlayers - 5) >> B.at(nlayers - 5) >> C.at(nlayers - 5) >> vs.at( 00083 nlayers - 5) >> E.at(nlayers - 5) >> densities.at(nlayers - 5); 00084 thicknesses.at(nlayers - 5) /= 1000.0; 00085 vp.at(nlayers - 5) /= 1000.0; 00086 vs.at(nlayers - 5) /= 1000.0; 00087 densities.at(nlayers - 5) /= 1000.0; 00088 infile.getline(title, titlelength); 00089 //for the layer n-3 00090 infile >> theta.at(nlayers - 5) >> phi.at(nlayers - 5); 00091 infile.getline(title, titlelength); 00092 infile >> thicknesses.at(nlayers - 5) >> vp.at(nlayers - 5) >> B.at(nlayers - 5) >> C.at(nlayers - 5) >> vs.at( 00093 nlayers - 5) >> E.at(nlayers - 5) >> densities.at(nlayers - 5); 00094 thicknesses.at(nlayers - 5) /= 1000.0; 00095 vp.at(nlayers - 5) /= 1000.0; 00096 vs.at(nlayers - 5) /= 1000.0; 00097 densities.at(nlayers - 5) /= 1000.0; 00098 infile.getline(title, titlelength); 00099 //for the layer n-2 00100 infile >> theta.at(nlayers - 5) >> phi.at(nlayers - 5); 00101 infile.getline(title, titlelength); 00102 infile >> thicknesses.at(nlayers - 5) >> vp.at(nlayers - 5) >> B.at(nlayers - 5) >> C.at(nlayers - 5) >> vs.at( 00103 nlayers - 5) >> E.at(nlayers - 5) >> densities.at(nlayers - 5); 00104 thicknesses.at(nlayers - 5) /= 1000.0; 00105 vp.at(nlayers - 5) /= 1000.0; 00106 vs.at(nlayers - 5) /= 1000.0; 00107 densities.at(nlayers - 5) /= 1000.0; 00108 infile.getline(title, titlelength); 00109 //for the layer n-1 : the last one is an half-space 00110 infile >> theta.at(nlayers - 5) >> phi.at(nlayers - 5); 00111 infile.getline(title, titlelength); 00112 infile >> thicknesses.at(nlayers - 5) >> vp.at(nlayers - 5) >> B.at(nlayers - 5) >> C.at(nlayers - 5) >> vs.at( 00113 nlayers - 5) >> E.at(nlayers - 5) >> densities.at(nlayers - 5); 00114 thicknesses.at(nlayers - 5) /= 1000.0; 00115 vp.at(nlayers - 5) /= 1000.0; 00116 vs.at(nlayers - 5) /= 1000.0; 00117 densities.at(nlayers - 5) /= 1000.0; 00118 infile.getline(title, titlelength); 00119 //the file contains depths we need thicknesses 00120 std::adjacent_difference(thicknesses.begin(), thicknesses.end(), 00121 thicknesses.begin()); 00122 } 00123 */ 00124 00125 // Function Read Model if uniform half-space 00126 /*void AnisoSurfaceWaveModel::ReadModel(const std::string &filename) 00127 { 00128 std::ifstream infile(filename.c_str()); 00129 size_t nlayers = 0; 00130 const size_t titlelength = 512; 00131 char title[titlelength]; 00132 //the first line of the model is the title that we don't care about 00133 infile.getline(title, titlelength); 00134 //the number of layers in the model does not include the bottom half-space 00135 infile >> nlayers; 00136 //so for us it is one layer more 00137 nlayers += 1; 00138 //there can be comments at the end of each line 00139 //so read the rest of the line, if nothing is there, it only reads the newline 00140 infile.getline(title, titlelength); 00141 //reserve memory for all the parameters 00142 // for a uniform in half-space 00143 thicknesses.assign(nlayers, 0.0); 00144 vp.assign(nlayers, 0.0); 00145 vs.assign(nlayers, 0.0); 00146 B.assign(nlayers, 0.0); 00147 C.assign(nlayers, 0.0); 00148 E.assign(nlayers, 0.0); 00149 theta.assign(nlayers, 0.0); 00150 phi.assign(nlayers, 0.0); 00151 densities.assign(nlayers, 0.0); 00152 for (size_t i = 0; i < nlayers; ++i) 00153 { 00154 //first line for each layer are the anisotropy angles 00155 infile >> theta.at(i) >> phi.at(i); 00156 //plus a comment 00157 infile.getline(title, titlelength); 00158 //then the other layer parameters 00159 infile >> thicknesses.at(i) >> vp.at(i) >> B.at(i) >> C.at(i) >> vs.at( 00160 i) >> E.at(i) >> densities.at(i); 00161 //transform from m to km for thickness and velocities and kg/m^3 to g/cm^3 for density 00162 thicknesses.at(i) /= 1000.0; 00163 vp.at(i) /= 1000.0; 00164 vs.at(i) /= 1000.0; 00165 densities.at(i) /= 1000.0; 00166 //swallow comment 00167 infile.getline(title, titlelength); 00168 } 00169 } 00170 */ 00171 00172 // Function Read Model if gradient in HS and Vs fixed at 4.87km/s at 410 km depth 00173 void AnisoSurfaceWaveModel::ReadModel(const std::string &filename) 00174 { 00175 std::ifstream infile(filename.c_str()); 00176 size_t nlayers = 0; 00177 int nhs = 5; 00178 const size_t titlelength = 512; 00179 char title[titlelength]; 00180 //the first line of the model is the title that we don't care about 00181 infile.getline(title, titlelength); 00182 //the number of layers in the model does not include the bottom half-space 00183 infile >> nlayers; 00184 //so for us it is one layer more 00185 nlayers += 1; 00186 //there can be comments at the end of each line 00187 //so read the rest of the line, if nothing is there, it only reads the newline 00188 infile.getline(title, titlelength); 00189 //reserve memory for all the parameters 00190 // for a uniform in half-space 00191 thicknesses.assign(nlayers, 0.0); 00192 vp.assign(nlayers, 0.0); 00193 vs.assign(nlayers, 0.0); 00194 B.assign(nlayers, 0.0); 00195 C.assign(nlayers, 0.0); 00196 E.assign(nlayers, 0.0); 00197 theta.assign(nlayers, 0.0); 00198 phi.assign(nlayers, 0.0); 00199 densities.assign(nlayers, 0.0); 00200 for (size_t i = 0; i < nlayers -1 ; ++i) 00201 { 00202 //first line for each layer are the anisotropy angles 00203 infile >> theta.at(i) >> phi.at(i); 00204 //plus a comment 00205 infile.getline(title, titlelength); 00206 //then the other layer parameters 00207 infile >> thicknesses.at(i) >> vp.at(i) >> B.at(i) >> C.at(i) >> vs.at( 00208 i) >> E.at(i) >> densities.at(i); 00209 //transform from m to km for thickness and velocities and kg/m^3 to g/cm^3 for density 00210 thicknesses.at(i) /= 1000.0; 00211 vp.at(i) /= 1000.0; 00212 vs.at(i) /= 1000.0; 00213 densities.at(i) /= 1000.0; 00214 //swallow comment 00215 infile.getline(title, titlelength); 00216 } 00217 } 00218 00219 // Function Write Model for a "gradient" in half-space 00220 00221 /*void AnisoSurfaceWaveModel::WriteModel(const std::string &filename) const 00222 { 00223 const size_t nlayers = thicknesses.size() +4; 00224 //check that all the layers have all parameters 00225 assert(nlayers> 0); 00226 assert(thicknesses.size() == vp.size()); 00227 assert(thicknesses.size() == vs.size()); 00228 assert(thicknesses.size() == B.size()); 00229 assert(thicknesses.size() == C.size()); 00230 assert(thicknesses.size() == E.size()); 00231 assert(thicknesses.size() == theta.size()); 00232 assert(thicknesses.size() == phi.size()); 00233 assert(thicknesses.size() == densities.size()); 00234 00235 std::ofstream outfile(filename.c_str()); 00236 //we have to write some title 00237 outfile << "Surface Wave model generated by gplib" << std::endl; 00238 //the code needs the number of layers without the half-space, we count the half-space 00239 outfile << nlayers -1 << std::endl; 00240 //we have to transform thicknesses to depths 00241 double currdepth = 0.0; 00242 for (size_t i = 0; i < nlayers -5 ; ++i) 00243 { 00244 //write the anisotropy angles, the new line is important because the original file format can have comments 00245 outfile << theta.at(i) << " " << phi.at(i) << "\n"; 00246 //calculate the depth to the bottom in m 00247 currdepth += thicknesses.at(i) * 1000.0; 00248 //write out the other parameters, again the newline matters 00249 //seismic anisotropic coefficient have to be positive (absolute value) 00250 outfile << currdepth << " " << vp.at(i) * 1000.0 << " " 00251 << fabs(B.at(i)) << " " << C.at(i) << " " << vs.at(i) * 1000.0 00252 << " " << fabs(E.at(i)) << " " << densities.at(i) * 1000.0 << "\n"; 00253 } 00254 // Divide the last layer (250 km) in 5 equal layers (50 km thick each) ( z>200 km) 00255 // we impose a gradient for Vs, Vp and densities 00256 outfile << theta.at(nlayers -5) << " " << phi.at(nlayers -5) << "\n"; 00257 currdepth += (thicknesses.at(nlayers -5)/5) * 1000.0; 00258 outfile << currdepth << " " << vp.at(nlayers -5) * 1000.0 << " " 00259 << fabs(B.at(nlayers -5)) << " " << C.at(nlayers -5) << " " << vs.at(nlayers -5) * 1000.0 00260 << " " << fabs(E.at(nlayers -5)) << " " << densities.at(nlayers -5) * 1000.0 << "\n"; 00261 // Multiply Vp by 1.02, Vs by 1.02 and densities by 1.01 00262 outfile << theta.at(nlayers -5) << " " << phi.at(nlayers -5) << "\n"; 00263 currdepth += (thicknesses.at(nlayers -5)/5) * 1000.0; 00264 outfile << currdepth << " " << vp.at(nlayers -5)*1.02 * 1000.0 << " " 00265 << fabs(B.at(nlayers -5)) << " " << C.at(nlayers -5) << " " << vs.at(nlayers -5)*1.02 * 1000.0 00266 << " " << fabs(E.at(nlayers -5)) << " " << densities.at(nlayers -5)*1.01 * 1000.0 << "\n"; 00267 // Mulpiply Vp by 1.04, Vs by 1.04, densities by 1.03 00268 outfile << theta.at(nlayers -5) << " " << phi.at(nlayers -5) << "\n"; 00269 currdepth += (thicknesses.at(nlayers -5)/5) * 1000.0; 00270 outfile << currdepth << " " << vp.at(nlayers -5)*1.04 * 1000.0 << " " 00271 << fabs(B.at(nlayers -5)) << " " << C.at(nlayers -5) << " " << vs.at(nlayers -5)*1.04 * 1000.0 00272 << " " << fabs(E.at(nlayers -5)) << " " << densities.at(nlayers -5)*1.03 * 1000.0 << "\n"; 00273 // Multiply Vp by 1.06, Vs by 1.06 and densities by 1.05 00274 outfile << theta.at(nlayers -5) << " " << phi.at(nlayers -5) << "\n"; 00275 currdepth += (thicknesses.at(nlayers -5)/5) * 1000.0; 00276 outfile << currdepth << " " << vp.at(nlayers -5)*1.06 * 1000.0 << " " 00277 << fabs(B.at(nlayers -5)) << " " << C.at(nlayers -5) << " " << vs.at(nlayers -5)*1.06 * 1000.0 00278 << " " << fabs(E.at(nlayers -5)) << " " << densities.at(nlayers -5)*1.05 * 1000.0 << "\n"; 00279 // Multiply Vp by 1.08, Vs by 1.08 and densities by 1.07 00280 outfile << theta.at(nlayers -5) << " " << phi.at(nlayers -5) << "\n"; 00281 currdepth += (thicknesses.at(nlayers -5)/5) * 1000.0; 00282 outfile << currdepth << " " << vp.at(nlayers -5)*1.08 * 1000.0 << " " 00283 << fabs(B.at(nlayers -5)) << " " << C.at(nlayers -5) << " " << vs.at(nlayers -5)*1.08 * 1000.0 00284 << " " << fabs(E.at(nlayers -5)) << " " << densities.at(nlayers -5)*1.07 * 1000.0 << "\n"; 00285 } 00286 */ 00287 00288 // Function Write Model if no gradient in half space 00289 /*void AnisoSurfaceWaveModel::WriteModel(const std::string &filename) const 00290 { 00291 const size_t nlayers = thicknesses.size(); 00292 //check that all the layers have all parameters 00293 assert(nlayers> 0); 00294 assert(thicknesses.size() == nlayers); 00295 assert(thicknesses.size() == nlayers); 00296 assert(thicknesses.size() == nlayers); 00297 assert(thicknesses.size() == nlayers); 00298 assert(thicknesses.size() == nlayers); 00299 assert(thicknesses.size() == nlayers); 00300 assert(thicknesses.size() == nlayers); 00301 assert(thicknesses.size() == nlayers); 00302 00303 std::ofstream outfile(filename.c_str()); 00304 //we have to write some title 00305 outfile << "Surface Wave model generated by gplib" << std::endl; 00306 //the code needs the number of layers without the half-space, we count the half-space 00307 outfile << nlayers -1 << std::endl; 00308 //we have to transform thicknesses to depths 00309 double currdepth = 0.0; 00310 for (size_t i = 0; i < nlayers; ++i) 00311 { 00312 //write the anisotropy angles, the new line is important because the original file format can have comments 00313 outfile << theta.at(i) << " " << phi.at(i) << "\n"; 00314 //calculate the depth to the bottom in m 00315 currdepth += thicknesses.at(i) * 1000.0; 00316 //write out the other parameters, again the newline matters 00317 //seismic anisotropic coefficient have to be positive (absolute value) 00318 outfile << currdepth << " " << vp.at(i) * 1000.0 << " " 00319 << fabs(B.at(i)) << " " << C.at(i) << " " << vs.at(i) * 1000.0 00320 << " " << fabs(E.at(i)) << " " << densities.at(i) * 1000.0 << "\n"; 00321 } 00322 } 00323 */ 00324 00325 // Function Write Model if gradient in HS and Vs fixed at 4.87km/s at 410 km depth 00326 void AnisoSurfaceWaveModel::WriteModel(const std::string &filename) const 00327 { 00328 int nhs = 5; // number of layers in half space 00329 const size_t nlayers = thicknesses.size(); 00330 //check that all the layers have all parameters 00331 assert(nlayers> 0); 00332 assert(thicknesses.size() == nlayers); 00333 assert(thicknesses.size() == nlayers); 00334 assert(thicknesses.size() == nlayers); 00335 assert(thicknesses.size() == nlayers); 00336 assert(thicknesses.size() == nlayers); 00337 assert(thicknesses.size() == nlayers); 00338 assert(thicknesses.size() == nlayers); 00339 assert(thicknesses.size() == nlayers); 00340 const double maxdepth = 410.; 00341 std::ofstream outfile(filename.c_str()); 00342 //we have to write some title 00343 outfile << "Surface Wave model generated by gplib" << std::endl; 00344 //the code needs the number of layers without the half-space, we count the half-space 00345 outfile << nlayers +nhs -2 << std::endl; 00346 //we have to transform thicknesses to depths 00347 double currdepth = 0.0; 00348 double thick = 0.0; 00349 double gradP = 0.0; 00350 double gradS = 0.0; 00351 double gradD = 0.0; 00352 int i = 0; 00353 for (size_t i = 0; i < nlayers -1; ++i) 00354 { 00355 //write the anisotropy angles, the new line is important because the original file format can have comments 00356 outfile << theta.at(i) << " " << phi.at(i) << "\n"; 00357 //calculate the depth to the bottom in m 00358 currdepth += thicknesses.at(i) * 1000.0; 00359 //write out the other parameters, again the newline matters 00360 //seismic anisotropic coefficient have to be positive (absolute value) 00361 outfile << currdepth << " " << vp.at(i) * 1000.0 << " " 00362 << fabs(B.at(i)) << " " << C.at(i) << " " << vs.at(i) * 1000.0 00363 << " " << fabs(E.at(i)) << " " << densities.at(i) * 1000.0 << "\n"; 00364 } 00365 thick = (maxdepth*1000 - currdepth)/nhs; 00366 gradS = ((vs.at(nlayers-1) - vs.at(nlayers-2))*1000)/(thick*nhs); 00367 gradP = ((vp.at(nlayers-1) - vp.at(nlayers-2))*1000)/(thick*nhs); 00368 gradD = ((densities.at(nlayers-1) - densities.at(nlayers-2))*1000)/(thick*nhs); 00369 for (i = 1; i < nhs+1 ; ++i) 00370 { 00371 outfile << theta.at(nlayers-1) << " " << phi.at(nlayers-1) << "\n"; 00372 currdepth += thick; 00373 outfile << currdepth << " " << vp.at(nlayers-2) * 1000.0 + gradP*thick*i << " " 00374 << fabs(B.at(nlayers-1)) << " " << C.at(nlayers-1) << " " << vs.at(nlayers-2) * 1000.0 + gradS*thick*i 00375 << " " << fabs(E.at(nlayers-1)) << " " << densities.at(nlayers-2) * 1000.0 + gradD*thick*i << "\n"; 00376 } 00377 } 00378 00379 void AnisoSurfaceWaveModel::WriteRunFile(const std::string &filename) const 00380 { 00381 std::ofstream runfile(filename.c_str()); 00382 boost::filesystem::path Path(filename); 00383 std::string dirname(filename + "_dir"); 00384 runfile << "#!/bin/bash" << std::endl; 00385 runfile << "mkdir " << dirname << std::endl; 00386 runfile << "mv " << filename << ".* " << dirname << std::endl; 00387 runfile << "cd " << dirname << std::endl; 00388 runfile << "cp " << Path.leaf() << ".mod animodel" << std::endl; 00389 // Run on Wegener 00390 runfile << "aniprop 2>&1> /dev/null" << std::endl; 00391 // Run on Leda or Stokes 00392 // runfile << "../bin/aniprop 2>&1> /dev/null" << std::endl; 00393 runfile << "cp out_cvel_R ../" << filename + ".cvel" << std::endl; 00394 if (!runfile.good()) 00395 throw FatalException("Problem writing to file: " + filename); 00396 } 00397 00398 // function Write Plot for a "gradient" in half-space 00399 00400 /*void AnisoSurfaceWaveModel::WritePlot(const std::string &filename) const 00401 { 00402 std::ofstream outfile(filename.c_str()); 00403 double cumthick = 0; 00404 const size_t nlayers = thicknesses.size() + 4; 00405 for (unsigned int i = 0; i < nlayers - 5 ; ++i) 00406 { 00407 //we write S-wave velocity and one anisotropy parameter and strike 00408 //once for the top 00409 outfile << cumthick << " " << vs.at(i) << " " << fabs(B.at(i)) << " " 00410 << phi.at(i) << std::endl; 00411 cumthick += thicknesses.at(i); 00412 //and once for the bottom 00413 outfile << cumthick << " " << vs.at(i) << " " << fabs(B.at(i)) << " " 00414 << phi.at(i) << std::endl; 00415 } 00416 // and then for the half-space divided in 5 layers 00417 // the first layer of the half-space (Vs back) 00418 //once for the top 00419 outfile << cumthick << " " << vs.back() << " " << fabs(B.back()) << " " 00420 << phi.back() << std::endl; 00421 cumthick += thicknesses.at(nlayers -5)/5; 00422 //and once for the bottom 00423 outfile << cumthick << " " << vs.back() << " " << fabs(B.back()) << " " 00424 << phi.back() << std::endl; 00425 // the second layer of the half-space (Vs back * 1.02) 00426 //once for the top 00427 outfile << cumthick << " " << vs.back()*1.02 << " " << fabs(B.back()) << " " 00428 << phi.back() << std::endl; 00429 cumthick += thicknesses.at(nlayers -5)/5; 00430 //and once for the bottom 00431 outfile << cumthick << " " << vs.back()*1.02 << " " << fabs(B.back()) << " " 00432 << phi.back() << std::endl; 00433 // the third layer of the half-space (Vs back * 1.04) 00434 //once for the top 00435 outfile << cumthick << " " << vs.back()*1.04 << " " << fabs(B.back()) << " " 00436 << phi.back() << std::endl; 00437 cumthick += thicknesses.at(nlayers -5)/5; 00438 //and once for the bottom 00439 outfile << cumthick << " " << vs.back()*1.04 << " " << fabs(B.back()) << " " 00440 << phi.back() << std::endl; 00441 // the fourth layer of the half-space (Vs back * 1.06) 00442 //once for the top 00443 outfile << cumthick << " " << vs.back()*1.06 << " " << fabs(B.back()) << " " 00444 << phi.back() << std::endl; 00445 cumthick += thicknesses.at(nlayers -5)/5; 00446 //and once for the bottom 00447 outfile << cumthick << " " << vs.back()*1.06 << " " << fabs(B.back()) << " " 00448 << phi.back() << std::endl; 00449 // the fifth layer of the half-space (Vs back * 1.08) 00450 //make sure we get a thick layer for the bottom half-space for plotting 00451 outfile << cumthick * 2.0 << " " << vs.back()*1.08 << " " << fabs(B.back()) << " " 00452 << phi.back() << std::endl; 00453 } 00454 */ 00455 00456 // Function Write Plot for a uniform half-space 00457 /*void AnisoSurfaceWaveModel::WritePlot(const std::string &filename) const 00458 { 00459 std::ofstream outfile(filename.c_str()); 00460 double cumthick = 0; 00461 const size_t nlayers = thicknesses.size(); 00462 for (unsigned int i = 0; i < nlayers ; ++i) 00463 { 00464 //we write S-wave velocity and one anisotropy parameter and strike 00465 //once for the top 00466 outfile << cumthick << " " << vs.at(i) << " " << fabs(B.at(i)) << " " 00467 << phi.at(i) << std::endl; 00468 cumthick += thicknesses.at(i); 00469 //and once for the bottom 00470 outfile << cumthick << " " << vs.at(i) << " " << fabs(B.at(i)) << " " 00471 << phi.at(i) << std::endl; 00472 } 00473 //make sure we get a thick layer for the bottom half-space for plotting 00474 outfile << cumthick * 2.0 << " " << vs.back() << " " << fabs(B.back()) << " " 00475 << phi.back() << std::endl; 00476 } 00477 */ 00478 00479 // Function Write Plot with gradient in HS and Vs fixed at 4.87km/s at 410km depth 00480 void AnisoSurfaceWaveModel::WritePlot(const std::string &filename) const 00481 { 00482 std::ofstream outfile(filename.c_str()); 00483 double cumthick = 0; 00484 int nhs = 5; // number of layers in half space 00485 const size_t nlayers = thicknesses.size(); 00486 for (unsigned int i = 0; i < nlayers -1 ; ++i) 00487 { 00488 //we write S-wave velocity and one anisotropy parameter and strike 00489 //once for the top 00490 outfile << cumthick << " " << vs.at(i) << " " << fabs(B.at(i)) << " " 00491 << phi.at(i) << std::endl; 00492 cumthick += thicknesses.at(i); 00493 //and once for the bottom 00494 outfile << cumthick << " " << vs.at(i) << " " << fabs(B.at(i)) << " " 00495 << phi.at(i) << std::endl; 00496 } 00497 const double maxdepth = 410.; 00498 double thick = 0.0; 00499 double gradS = 0.0; 00500 int i = 0; 00501 thick = (maxdepth - cumthick)/nhs; 00502 gradS = (vs.at(nlayers-1) - vs.at(nlayers-2))/(thick*nhs); 00503 // the first layer of the half-space (Vs back) 00504 for (i = 1; i < nhs +1; ++i) 00505 { 00506 //once for the top 00507 outfile << cumthick << " " << vs.at(nlayers-2) + gradS*thick*i << " " 00508 << fabs(B.at(nlayers-1)) << " " << phi.at(nlayers-1)<< std::endl; 00509 cumthick += thick; 00510 outfile << cumthick << " " << vs.at(nlayers-2) + gradS*thick*i << " " 00511 << fabs(B.at(nlayers-1)) << " " << phi.at(nlayers-1)<< std::endl; 00512 } 00513 //make sure we get a thick layer for the bottom half-space for plotting 00514 outfile << cumthick * 2.0 << " " << vs.back() << " " << fabs(B.back()) << " " 00515 << phi.back() << std::endl; 00516 } 00517 }
1.5.8