GPLIB++
FkModel.cpp
Go to the documentation of this file.
1 #include "FkModel.h"
2 
3 #include <fstream>
4 #include <iomanip>
5 using namespace std;
6 namespace gplib
7  {
8  FkModel::FkModel():
10  {
11  }
13  {
14  }
15 
16  void FkModel::WriteModel(const std::string filename)
17  {
18  ofstream neu(filename.c_str()); //the file to write the model to
19  int sourceindex = 0; // the index of the source layer
20  double sourcethickness = 0; //thickness of the layer above the source (FK needs dummy layer)
21  double totalthickness = 0; // counter to add up thicknesses to calculate sourcethickness
22  unsigned int totallayers = GetThickness().size();
23 
24  sourceindex = FindLayer(GetSourceDepth()); //find out in which layer the source lies
25  for (int i = 0; i < sourceindex; ++i)
26  totalthickness += GetThickness().at(i); // add up the thickness of all layers above the source (without sourcelayer)
27  sourcethickness = GetSourceDepth() - totalthickness; //calculate the thickness of the sourcelayer so that source lies
28  // at lower boundary
29  neu
30  << ".F. \n"; //write File header
31  neu
32  << " 0 64 \n"; //FK-Code is VERY
33  neu << "GREEN.1\n"; //picky about spaces etc. !!
34  neu << " 6.0 " << setprecision(1) << setiosflags(
35  ios::scientific) << GetSourceDepth() << resetiosflags(
36  ios::scientific) << " 1" << setw(5) << setfill(' ')
37  << setprecision(8) << GetNpts() / 2 << setw(5) << setfill(' ')
38  << GetNpts() << " " << setprecision(3)
39  << setiosflags(ios::fixed) << GetDt() << resetiosflags(ios::fixed)
40  << setfill(' ') << setw(6) << totallayers << " 1\n"; //see documentation
41  // of FK for meaning
42  neu
43  << " 1 1 1 1 1 1 1 1 1 1 0 \n";
44 
45  for (unsigned int i = 0; i < totallayers; ++i) //write all layers
46  {
47  neu << setprecision(4) << setiosflags(ios::scientific) << " "
48  << GetThickness().at(i) << " " << GetPVelocity().at(i) << " ";
49  neu << GetSVelocity().at(i) << " " << GetDensity().at(i) << " "
50  << GetQp().at(i) << " " << GetQs().at(i) << endl;
51  }
52 
53  neu << " " << sourceindex + 1
54  << " \n";
55  neu << " 0.4000000E+03 1.500000E+00 0\n"; //write epilogue for FK-Code
56  neu << setw(5) << setfill(' ') << GetRecDist().size()
57  << " 10000.0 30.0 0.4 0.3\n";
58  for (unsigned int i = 0; i < GetRecDist().size(); ++i)
59  neu << " " << setw(3) << setprecision(3) << setiosflags(ios::fixed
60  | ios::showpoint) << GetRecDist().at(i)
61  << " 0.0 9999.0 \n";
62  neu << " " << flush;
63  neu << "\n";
64  }
65 
66  void FkModel::ReadModel(const std::string filename)
67  {
68  ifstream infile(filename.c_str());
69  char dummy[255];
70  double number;
71  int nlayers, nstat;
72  double SourceDepth, dt;
73  int npts;
74  infile.getline(dummy, 255);
75  infile.getline(dummy, 255);
76  infile.getline(dummy, 255);
77  infile >> number >> SourceDepth >> number >> number >> npts >> dt
78  >> nlayers >> number;
79  infile.getline(dummy, 255);
80  infile.getline(dummy, 255);
81  SetSourceDepth(SourceDepth);
82  SetNpts(npts);
83  SetDt(dt);
84  Init(nlayers);
85  for (int i = 0; i < nlayers; ++i)
86  infile >> SetThickness().at(i) >> SetPVelocity().at(i)
87  >> SetSVelocity().at(i) >> SetDensity().at(i) >> SetQp().at(i)
88  >> SetQs().at(i);
89  infile >> number;
90  infile.getline(dummy, 255);
91  infile.getline(dummy, 255);
92  infile >> nstat >> number >> number >> number >> number;
93  SetRecDist().assign(nstat, 0);
94  for (int i = 0; i < nstat; ++i)
95  infile >> SetRecDist().at(i) >> number >> number;
96  }
97  }
int FindLayer(const double depth)
Returns the index of the layer that correponds to depth in km.
trealdata & SetDensity()
Read-write access to the vector of densities in g/cm^3 in each layer.
Definition: SeismicModel.h:96
double GetSourceDepth() const
Get the depth to the seismic source that generates the wavefield.
Definition: SeismicModel.h:51
virtual void ReadModel(const std::string filename)
Read the model in its native format from a file.
Definition: FkModel.cpp:66
const trealdata & GetDensity() const
Read-only access to the vector of densities in g/cm^3 in each layer.
Definition: SeismicModel.h:91
trealdata & SetQs()
Read-write access to the vector of S quality factors for each layer.
Definition: SeismicModel.h:126
virtual void WriteModel(const std::string filename)
Write the model in its native format to a file.
Definition: FkModel.cpp:16
const trealdata & GetQp() const
Read-only access to the vector of P quality factors for each layer.
Definition: SeismicModel.h:111
trealdata & SetThickness()
Read-write access to the vector of thicknesses in km in each layer.
Definition: SeismicModel.h:106
trealdata & SetQp()
Read-write access to the vector of P quality factors for each layer.
Definition: SeismicModel.h:116
void SetNpts(const unsigned int s)
Set the number of points for synthetic seismogram calculation.
Definition: SeismicModel.h:46
trealdata & SetPVelocity()
Read-write access to the vector of P-velocities in km/s in each layer.
Definition: SeismicModel.h:76
The class SeismicModel is the base class for some of the model format for seismic codes...
Definition: SeismicModel.h:17
virtual ~FkModel()
Definition: FkModel.cpp:12
void SetSourceDepth(const double s)
Set the depth to the seismic source.
Definition: SeismicModel.h:56
double GetDt() const
Get the time between two samples in s, this is for synthetic forward calculations.
Definition: SeismicModel.h:61
const trealdata & GetSVelocity() const
Read-only access to the vector of S-velocities in km/s in each layer.
Definition: SeismicModel.h:81
void Init(const int nlayers)
Init provides a convenient way to allocate memory in all structures for a given number of layers...
trealdata & SetRecDist()
Read-write access to the vector of receiver distances.
Definition: FkModel.h:23
trealdata & SetSVelocity()
Read-write access to the vector of S-velocities in km/s in each layer.
Definition: SeismicModel.h:86
unsigned int GetNpts() const
Get the number of points for synthetic seismogram calculation.
Definition: SeismicModel.h:41
const trealdata & GetRecDist() const
Read-only access to the vector of receiver distances.
Definition: FkModel.h:18
void SetDt(const double s)
Set the time between two samples in s, this is for synthetic forward calculations.
Definition: SeismicModel.h:66
const trealdata & GetPVelocity() const
Read-only access to the vector of P-velocities in km/s in each layer.
Definition: SeismicModel.h:71
const trealdata & GetThickness() const
Read-only access to the vector of thicknesses in km in each layer.
Definition: SeismicModel.h:101
const trealdata & GetQs() const
Read-only access to the vector of S quality factors for each layer.
Definition: SeismicModel.h:121