GPLIB++
AnisoSurfaceWaveModel.cpp
Go to the documentation of this file.
1 #include "AnisoSurfaceWaveModel.h"
2 #include <fstream>
3 #include <cassert>
4 #include <algorithm>
5 #include <numeric>
6 #include <boost/bind.hpp>
7 #include <string>
8 #include <boost/filesystem.hpp>
9 #include "FatalException.h"
10 #include <iostream>
11 
12 namespace gplib
13  {
14 
15  AnisoSurfaceWaveModel::AnisoSurfaceWaveModel()
16  {
17  }
18 
19  AnisoSurfaceWaveModel::~AnisoSurfaceWaveModel()
20  {
21  }
22 
23  void AnisoSurfaceWaveModel::ReadModel(const std::string &filename)
24  {
25  std::cout << "Read Model " << std::endl;
26  }
27 
28  // Function Read Model for MODES
29  /*void AnisoSurfaceWaveModel::ReadModel(const std::string &filename)
30  {
31  std::ifstream infile(filename.c_str());
32  size_t nlayers = 0;
33  const size_t titlelength = 512;
34  char title[titlelength];
35  //the first line of the model is the title that we don't care about
36  infile.getline(title, titlelength);
37  //the second layer is fixed parameters we don't care about
38  infile >> 0 >> 50 >> 1;
39  //the number of layers in the model does include the bottom fixed model (14 layers) + layers for boundaries
40  infile >> nlayers >> 3 >> 5;
41  //so the number of inverted layers
42  nlayers = (nlayers - 15)/2 + 1;
43  //there can be comments at the end of each line
44  //so read the rest of the line, if nothing is there, it only reads the newline
45  infile.getline(title, titlelength);
46  //reserve memory for all the parameters
47  // for a uniform in half-space
48  thicknesses.assign(nlayers, 0.0);
49  vp.assign(nlayers, 0.0);
50  vsapp.assign(nlayers, 0.0);
51  B.assign(nlayers, 0.0);
52  C.assign(nlayers, 0.0);
53  E.assign(nlayers, 0.0);
54  theta.assign(nlayers, 0.0);
55  phi.assign(nlayers, 0.0);
56  densities.assign(nlayers, 0.0);
57  //we read the fixed bottom part of the model
58  infile >> 0 >> 13012.20 >> 11218.76 >> 3613.67 >> 99999 >> 85 >> 6371.00;
59  infile >> 101430 >> 13010.00 >> 11217.18 >> 3612.58 >> 99999 >> 85 >> 6269.57;
60  infile >> 1217500 >> 12703.70 >> 11002.11 >> 3452.58 >> 99999 >> 85 >> 5153.50;
61  infile >> 1217500 >> 12139.10 >> 10288.78 >> 0.00 >> 99999 >> 0.0 >> 5153.50;
62  infile >> 3479500 >> 9914.50 >> 7999.83 >> 0.00 >> 99999 >> 0.0 >> 2891.50;
63  infile >> 3479500 >> 5772.10 >> 13621.92 >> 7248.53 >> 99999 >> 274.0 >> 2891.50;
64  infile >> 3878500 >> 5331.30 >> 13330.80 >> 7112.86 >> 99999 >> 371.8 >> 2492.50;
65  infile >> 4373500 >> 5102.70 >> 12771.47 >> 6898.60 >> 99999 >> 414.9 >> 1997.50;
66  infile >> 4868500 >> 4856.20 >> 12170.18 >> 6663.31 >> 99999 >> 463.0 >> 1502.50;
67  infile >> 5363500 >> 4592.60 >> 11452.03 >> 6370.39 >> 99999 >> 510.2 >> 1007.50;
68  infile >> 5711000 >> 4238.70 >> 10775.04 >> 5947.18 >> 99999 >> 549.5 >> 660.00;
69  infile >> 5711000 >> 3920.10 >> 10152.58 >> 5569.85 >> 99999 >> 172.9 >> 660.00;
70  infile >> 5961000 >> 3931.70 >> 9314.71 >> 5041.52 >> 99999 >> 162.5 >> 410.00;
71  infile >> 5961000 >> 3467.98 >> 8435.09 >> 4870.00 >> 99999 >> 140.0 >> 410.00;
72  //then we read the boundary between layer n-1 and fixed part of the model
73  double currdepth =0.0;
74  double currdepth_center = 0.0;
75  infile >> currdepth_center >> densities.at(nlayers-1) >> vp.at(nlayers-1) >> vsapp.at(nlayers-1)
76  >> 99999 >> 100 >> currdepth;
77  thicknesses(nlayers-1) = 410.00 - currdepth;
78  //finally we read the inverted part of the model
79  int i = nlayers -2;
80  while (i >= 0)
81  {
82  infile >> currdepth_center >> densities.at(i) >> vp.at(i) >> vsapp.at(i) >> 99999
83  >> 100 >> currdepth;
84  infile >> currdepth_center + thicknesses.at(i) >> densities.at(i) >> vp.at(i) >> vsapp.at(i) >> 99999
85  >> 100 >> currdepth - thicknesses.at(i)/1000;
86  i--;
87  }
88  //transform from m to km for thickness and velocities and kg/m^3 to g/cm^3 for density
89  vp.at(i) /= 1000.0;
90  vsapp.at(i) /= 1000.0;
91  densities.at(i) /= 1000.0;
92  thicknesses.at(i) /=1000.0;
93  }
94  */
95 
96  // Function Write Model for MODES
97  void AnisoSurfaceWaveModel::WriteModel(const std::string &filename) const
98  {
99  const size_t nlayers = thicknesses.size();
100  //check that all the layers have all parameters
101  assert(nlayers> 0);
102  assert(thicknesses.size() == nlayers);
103  assert(vp.size() == nlayers);
104  assert(vsapp.size() == nlayers);
105  assert(densities.size() == nlayers);
106  assert(B.size() == nlayers);
107  assert(C.size() == nlayers);
108  assert(phi.size() == nlayers);
109  std::ofstream outfile(filename.c_str());
110  //we have to write some title
111  outfile << "title" << std::endl;
112  //the code needs some parameters (reference freq = 50 Hz)
113  outfile << 0 << " " << 50 << " " << 1 << " " << std::endl;
114  //the code needs the numbers of layers in the model / number of lay. in inner / outer core
115  outfile << (nlayers - 1) * 2 + 1 + 14 << " " << 3 << " " << 5 << " "
116  << std::endl;
117  //Model fixed from 410 km to the center of the Earth
118  outfile << 0 << " " << 13012.20 << " " << 11218.76 << " " << 3613.67
119  << " " << 99999 << " " << 85 << " " << 6371.00 << " " << "\n";
120  outfile << 101430 << " " << 13010.00 << " " << 11217.18 << " "
121  << 3612.58 << " " << 99999 << " " << 85 << " " << 6269.57 << " "
122  << "\n";
123  outfile << 1217500 << " " << 12703.70 << " " << 11002.11 << " "
124  << 3452.58 << " " << 99999 << " " << 85 << " " << 5153.50 << " "
125  << "\n";
126  outfile << 1217500 << " " << 12139.10 << " " << 10288.78 << " " << 0.00
127  << " " << 99999 << " " << 0.0 << " " << 5153.50 << " " << "\n";
128  outfile << 3479500 << " " << 9914.50 << " " << 7999.83 << " " << 0.00
129  << " " << 99999 << " " << 0.0 << " " << 2891.50 << " " << "\n";
130  outfile << 3479500 << " " << 5772.10 << " " << 13621.92 << " "
131  << 7248.53 << " " << 99999 << " " << 274.0 << " " << 2891.50 << " "
132  << "\n";
133  outfile << 3878500 << " " << 5331.30 << " " << 13330.80 << " "
134  << 7112.86 << " " << 99999 << " " << 371.8 << " " << 2492.50 << " "
135  << "\n";
136  outfile << 4373500 << " " << 5102.70 << " " << 12771.47 << " "
137  << 6898.60 << " " << 99999 << " " << 414.9 << " " << 1997.50 << " "
138  << "\n";
139  outfile << 4868500 << " " << 4856.20 << " " << 12170.18 << " "
140  << 6663.31 << " " << 99999 << " " << 463.0 << " " << 1502.50 << " "
141  << "\n";
142  outfile << 5363500 << " " << 4592.60 << " " << 11452.03 << " "
143  << 6370.39 << " " << 99999 << " " << 510.2 << " " << 1007.50 << " "
144  << "\n";
145  outfile << 5711000 << " " << 4238.70 << " " << 10775.04 << " "
146  << 5947.18 << " " << 99999 << " " << 549.5 << " " << 660.00 << " "
147  << "\n";
148  outfile << 5711000 << " " << 3920.10 << " " << 10152.58 << " "
149  << 5569.85 << " " << 99999 << " " << 172.9 << " " << 660.00 << " "
150  << "\n";
151  outfile << 5961000 << " " << 3931.70 << " " << 9314.71 << " "
152  << 5041.52 << " " << 99999 << " " << 162.5 << " " << 410.00 << " "
153  << "\n";
154  // we fix Vs=4.87 km/s at 410 km in Germany
155  outfile << 5961000 << " " << 3467.98 << " " << 8435.09 << " "
156  << 4870.00 << " " << 99999 << " " << 140.0 << " " << 410.00 << " "
157  << "\n";
158  // we fix Vs=5.00 km/s at 410 km depth in Slave
159  // outfile << 5961000 << " " << 3540.00 << " " << 9000.00 << " " << 5000.00 << " " << 99999 << " "
160  // << 140.0 << " " << 410.00 << " " << "\n";
161 
162  //we compute currdepth = maximum depth inverted
163  double currdepth = 0.0;
164  for (size_t j = 0; j < nlayers - 1; j++)
165  {
166  currdepth += thicknesses.at(j);
167  }
168  double currdepth_center = (6371 - currdepth) * 1000;
169  //the layer n-1 is not inverted for the seismic model: the thickness is fixed to 410 - currdepth
170  //and we impose a gradient between Vs(nlayer - 2) and 4.87km/s (fixed at 410 km depth)
171  outfile << currdepth_center << " " << (vs.at(nlayers - 2) * 0.554
172  + 0.77) * 1000 << " " << vs.at(nlayers - 2) * 1.80 * 1000 << " "
173  << vs.at(nlayers - 2) * 1000.0 << " " << 99999 << " " << 100 << " "
174  << currdepth << " " << "\n";
175 
176  int i = nlayers - 2;
177  while (i >= 0)
178  {
179  outfile << currdepth_center << " " << densities.at(i) * 1000.0
180  << " " << vp.at(i) * 1000 << " " << vsapp.at(i) * 1000.0 << " "
181  << 99999 << " " << 100 << " " << currdepth << " " << "\n";
182  //calculate the distance to the center of the Earth in m
183  currdepth -= thicknesses.at(i);
184  //calculate the distance to the surface in km
185  currdepth_center += thicknesses.at(i) * 1000;
186  //write out the different parameters
187  outfile << currdepth_center << " " << densities.at(i) * 1000.0
188  << " " << vp.at(i) * 1000 << " " << vsapp.at(i) * 1000.0 << " "
189  << 99999 << " " << 100 << " " << currdepth << " " << "\n";
190  i--;
191  }
192  }
193 
194  // Function WriteRunFile for Modes
195  void AnisoSurfaceWaveModel::WriteRunFile(const std::string &filename) const
196  {
197  std::ofstream runfile(filename.c_str());
198  std::ofstream inputfile((filename + ".in").c_str());
199  inputfile << "isomodel" << " " << std::endl;
200  inputfile << "output" << " " << std::endl;
201  inputfile << "phase" << " " << std::endl;
202  // periods for Germany and synthetics
203  inputfile << 34 << " " << std::endl;
204  inputfile << 14.08 << " " << 14.49 << " " << 14.92 << " " << 15.38
205  << " " << 15.87 << " " << 16.39 << " " << 16.94 << " " << 17.54
206  << " " << 18.18 << " " << 18.86 << " " << 19.60 << " " << 20.40
207  << " " << 21.27 << " " << 22.22 << " " << 23.25 << " " << 24.39
208  << " " << 25.64 << " " << 27.02 << " " << 29.41 << " " << 31.25
209  << " " << 33.33 << " " << 35.71 << " " << 38.46 << " " << 41.66
210  << " " << 45.45 << " " << 50. << " " << 55.55 << " " << 62.50
211  << " " << 71.42 << " " << 83.33 << " " << 100. << " " << 125.
212  << " " << 166.66 << " " << 250. << " " << std::endl;
213  // periods for Slave
214  /* inputfile << 31 << " " << std::endl;
215  inputfile << 15.15 << " " << 15.62 << " " << 16.12 << " " << 16.67 << " " << 17.24 << " " << 17.86
216  << " " << 18.52 << " " << 19.23 << " " << 20.00 << " " << 20.83 << " " << 21.74 << " " << 22.73
217  << " " << 23.81 << " " << 25.00 << " " << 26.31 << " " << 27.78 << " " << 29.41 << " " << 31.25
218  << " " << 33.33 << " " << 35.71 << " " << 38.46 << " " << 41.67 << " " << 45.45 << " " << 50.00
219  << " " << 55.55 << " " << 62.50 << " " << 71.40 << " " << 83.33 << " " << 100. << " " << 125.
220  << " " << 166.67 << " " << std::endl;
221  */
222  boost::filesystem::path Path(filename);
223  std::string dirname(filename + "_dir");
224  runfile << "#!/bin/bash" << std::endl;
225  runfile << "mkdir " << dirname << std::endl;
226  runfile << "mv " << filename << ".* " << dirname << std::endl;
227  runfile << "cd " << dirname << std::endl;
228  runfile << "cp " << Path.leaf() << ".mod isomodel" << std::endl;
229  runfile << "cp " << Path.leaf() << ".in input" << std::endl;
230  // Run on Wegener
231  runfile << "isoraylC_fund < input 2>&1> /dev/null" << std::endl;
232  runfile << "cp �MXAe ../" << filename + ".cvel" << std::endl;
233  // Run on Stokes
234  // runfile << "../bin/isoraylC_fund < input 2>&1> /dev/null" << std::endl;
235  // runfile << "cp phase ../" << filename + ".cvel" << std::endl;
236  if (!runfile.good())
237  throw FatalException("Problem writing to file: " + filename);
238  }
239 
240  //Function WritePlot for Modes
241  void AnisoSurfaceWaveModel::WritePlot(const std::string &filename) const
242  {
243  std::ofstream outfile(filename.c_str());
244  const size_t nlayers = thicknesses.size();
245  outfile << 3613.67 / 1000 << " " << 0 << " " << 0 << " " << 6371.00
246  << " " << "\n";
247  outfile << 3612.58 / 1000 << " " << 0 << " " << 0 << " " << 6269.57
248  << " " << "\n";
249  outfile << 3452.58 / 1000 << " " << 0 << " " << 0 << " " << 5153.50
250  << " " << "\n";
251  outfile << 0.00 << " " << 0 << " " << 0 << " " << 5153.50 << " "
252  << "\n";
253  outfile << 0.00 << " " << 0 << " " << 0 << " " << 2891.50 << " "
254  << "\n";
255  outfile << 7248.53 / 1000 << " " << 0 << " " << 0 << " " << 2891.50
256  << " " << "\n";
257  outfile << 7112.86 / 1000 << " " << 0 << " " << 0 << " " << 2492.50
258  << " " << "\n";
259  outfile << 6898.60 / 1000 << " " << 0 << " " << 0 << " " << 1997.50
260  << " " << "\n";
261  outfile << 6663.31 / 1000 << " " << 0 << " " << 0 << " " << 1502.50
262  << " " << "\n";
263  outfile << 6370.39 / 1000 << " " << 0 << " " << 0 << " " << 1007.50
264  << " " << "\n";
265  outfile << 5947.18 / 1000 << " " << 0 << " " << 0 << " " << 660.00
266  << " " << "\n";
267  outfile << 5569.85 / 1000 << " " << 0 << " " << 0 << " " << 660.00
268  << " " << "\n";
269  outfile << 5041.52 / 1000 << " " << 0 << " " << 0 << " " << 410.00
270  << " " << "\n";
271  outfile << 4870.00 / 1000 << " " << 0 << " " << 0 << " " << 410.00
272  << " " << "\n"; // Germany
273  // outfile << 5000.00/1000 << " " << 0 << " " << 0 << " " << 410.00 << " " << "\n"; // Slave
274  int j = 0;
275  //we compute currdepth = maximum depth inverted
276  double currdepth = 0.0;
277  for (j = 0; j < nlayers - 1; j++)
278  {
279  currdepth += thicknesses.at(j);
280  }
281  double currdepth_center = (6371 - currdepth) * 1000;
282  //one boundary for last inverted layer (thickness is not inverted but fixed to 410 - currdepth)
283  //we have a gradient between currdepth and 410 km depth
284  outfile << vs.at(nlayers - 2) << " " << 0 << " " << 0 << " "
285  << currdepth << " " << "\n";
286  int i = nlayers - 2;
287  while (i >= 0)
288  {
289  outfile << vs.at(i) << " " << B.at(i) << " " << phi.at(i) << " "
290  << currdepth << " " << "\n";
291  //calculate the distance to the center of the Earth in m
292  currdepth -= thicknesses.at(i);
293  //calculate the distance to the surface in km
294  currdepth_center += thicknesses.at(i) * 1000;
295  //write out the different parameters
296  outfile << vs.at(i) << " " << B.at(i) << " " << phi.at(i) << " "
297  << currdepth << " " << "\n";
298  i--;
299  }
300  }
301 
302  }
303 
gplib::AnisoSurfaceWaveModel::WriteModel
void WriteModel(const std::string &filename) const
Write the model to a file with name filename.
Definition: AnisoSurfaceWaveModel.cpp:97
gplib::AnisoSurfaceWaveModel::ReadModel
void ReadModel(const std::string &filename)
Read the model information from an ascii file with name filename.
Definition: AnisoSurfaceWaveModel.cpp:23
gplib::AnisoSurfaceWaveModel::WriteRunFile
void WriteRunFile(const std::string &filename) const
Write out a script that computes synthetic data for a model with name filename+.dat.
Definition: AnisoSurfaceWaveModel.cpp:195
gplib::AnisoSurfaceWaveModel::WritePlot
void WritePlot(const std::string &filename) const
Write out an ascii file for plotting the model with xmgrace.
Definition: AnisoSurfaceWaveModel.cpp:241
gplib::FatalException
The basic exception class for all errors that arise in gplib.
Definition: FatalException.h:11
Generated on Wed May 25 2016 14:21:06 for GPLIB++ by   doxygen 1.8.6