makeinput.cpp
Go to the documentation of this file.00001 #include <iostream>
00002 #include "C1DMTSynthData.h"
00003 #include "CResPkModel.h"
00004 #include "CRecCalc.h"
00005 #include "CUniformRNG.h"
00006 #include "MakeInputConf.h"
00007 #include "CSeismicDataComp.h"
00008 #include "Adaptors.h"
00009 #include <vector>
00010 #include <boost/bind.hpp>
00011 #include <string>
00012
00013 using namespace std;
00014 using namespace gplib;
00015
00016 string version = "$Id: makeinput.cpp 1816 2009-09-07 11:28:35Z mmoorkamp $";
00017
00018 int main(int argc, char *argv[])
00019 {
00020 cout
00021 << " This is makeinput: Create MT and receiver function models and data"
00022 << endl;
00023 cout << " to be used as input for inversion methods." << endl;
00024 cout
00025 << " If one argument is given, this argument is taken as the filename base for all output files."
00026 << endl;
00027 cout
00028 << " Otherwise the output file base is read from the configuration file 'makeinput.conf' "
00029 << endl;
00030 cout << " The configuration file contains all other relevant settings. "
00031 << endl;
00032 cout << " This is Version: " << version << endl << endl;
00033 UniformRNG Random;
00034 MakeInputConf Configuration;
00035 Configuration.GetData("makeinput.conf");
00036
00037 string outputbase;
00038 if (argc > 1)
00039 outputbase = argv[1];
00040 else
00041 outputbase = Configuration.outputbase;
00042
00043 C1DMTSynthData MTSynth;
00044 ResPkModel RecModel;
00045 RecCalc RecSynth(Configuration.shift, Configuration.omega,
00046 Configuration.sigma, Configuration.cc);
00047 SeismicDataComp RecSynthData;
00048 const double basevel = 3.2;
00049 const double basedens = 2.6;
00050 const double baseres = 1.0 + Random.GetNumber();
00051 const double basethick = 1.0 + Random.GetNumber() * 10;
00052 vector<double> velocities(Configuration.nlayers);
00053 vector<double> densities(Configuration.nlayers);
00054 vector<double> resistivities(Configuration.nlayers);
00055 vector<double> mtthickness(Configuration.nlayers);
00056 vector<double> recthickness(Configuration.nlayers);
00057 densities.front() = basedens;
00058 resistivities.front() = baseres;
00059 velocities.front() = basevel;
00060 mtthickness.front() = basethick;
00061 if (Configuration.correlated)
00062 {
00063 recthickness.front() = basethick;
00064 }
00065 else
00066 {
00067 recthickness.front() = basethick * Random.GetNumber(1.1, 2.0);
00068 }
00069 for (int i = 1; i < Configuration.nlayers; ++i)
00070 {
00071 mtthickness.at(i) = mtthickness.at(i - 1) * Random.GetNumber(1.1, 2.0);
00072 if (Configuration.correlated)
00073 {
00074 recthickness.at(i) = mtthickness.at(i);
00075 }
00076 else
00077 {
00078 recthickness.at(i) = mtthickness.at(i) * Random.GetNumber(1.1, 2.0);
00079
00080 }
00081 densities.at(i) = densities.at(i - 1) * Random.GetNumber(1.1, 1.5);
00082 resistivities.at(i) = Random.GetNumber(0.0, 3.0);
00083 velocities.at(i) = velocities.at(i - 1) * Random.GetNumber(1.05, 1.2);
00084 }
00085 transform(resistivities.begin(), resistivities.end(),
00086 resistivities.begin(), boost::bind(fopow<double, double> (), 10, _1));
00087 MTSynth.SetResistivities(resistivities);
00088 MTSynth.SetThicknesses(mtthickness);
00089 MTSynth.writemodel(outputbase + "_mt.mod");
00090 MTSynth.writeplot(outputbase + "_mt.plot");
00091 MTSynth.GetData();
00092 MTSynth.WriteAsMtt(outputbase);
00093
00094 RecModel.Init(Configuration.nlayers);
00095 RecModel.slowness = Configuration.slowness;
00096 copy(recthickness.begin(), recthickness.end(), RecModel.Thickness.begin());
00097 copy(velocities.begin(), velocities.end(), RecModel.SVelocity.begin());
00098 copy(densities.begin(), densities.end(), RecModel.Density.begin());
00099 transform(RecModel.SVelocity.begin(), RecModel.SVelocity.end(),
00100 RecModel.PVelocity.begin(), boost::bind(multiplies<double> (), _1,
00101 Configuration.poisson));
00102 RecModel.dt = Configuration.dt;
00103 RecModel.npts = Configuration.npts;
00104 RecSynth.CalcRecSynth(Configuration.outputbase + "_rec", RecModel,
00105 RecSynthData, false);
00106 RecModel.WritePlot(outputbase + "_rec");
00107 RecSynthData.WriteAsSac(outputbase + ".rec");
00108 }