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 string version = "$Id: makeinput.cpp 1093 2007-05-08 17:28:32Z max $";
00015 
00016 int main(int argc, char *argv[])
00017 {
00018         cout << " This is makeinput: Create MT and receiver function models and data" << endl; // write some info
00019         cout << " to be used as input for inversion methods." << endl;
00020         cout << " If one argument is given, this argument is taken as the filename base for all output files." << endl;
00021         cout << " Otherwise the output file base is read from the configuration file 'makeinput.conf' " << endl;
00022         cout << " The configuration file contains all other relevant settings. " << endl;
00023         cout << " This is Version: " << version << endl << endl;
00024         UniformRNG Random;
00025         MakeInputConf Configuration;
00026         Configuration.GetData("makeinput.conf");
00027         
00028         string outputbase;
00029         if (argc > 1)
00030                 outputbase = argv[1];
00031         else
00032                 outputbase = Configuration.outputbase;
00033                 
00034         C1DMTSynthData MTSynth;
00035         ResPkModel RecModel;
00036         RecCalc RecSynth(Configuration.shift, Configuration.omega, Configuration.sigma, Configuration.cc);
00037         SeismicDataComp RecSynthData;
00038         const double basevel = 3.2;
00039         const double basedens = 2.6;
00040         const double baseres = 1.0 + Random.GetNumber();
00041         const double basethick = 1.0 + Random.GetNumber()*10;
00042         vector<double> velocities(Configuration.nlayers);
00043         vector<double> densities(Configuration.nlayers);
00044         vector<double> resistivities(Configuration.nlayers);
00045         vector<double> mtthickness(Configuration.nlayers);
00046         vector<double> recthickness(Configuration.nlayers);
00047         densities.front() = basedens;
00048         resistivities.front() = baseres;
00049         velocities.front() = basevel;
00050         mtthickness.front() = basethick;
00051         if (Configuration.correlated)
00052         {
00053                 recthickness.front() = basethick;
00054         }
00055         else
00056         {
00057                 recthickness.front() = basethick * Random.GetNumber(1.1,2.0);
00058         }
00059         for (int i = 1; i < Configuration.nlayers; ++i)
00060         {
00061                 mtthickness.at(i) = mtthickness.at(i-1) * Random.GetNumber(1.1,2.0);
00062                 if (Configuration.correlated)
00063                 {
00064                          recthickness.at(i) = mtthickness.at(i);
00065                 }
00066                 else
00067                 {
00068                                 recthickness.at(i) = mtthickness.at(i) * Random.GetNumber(1.1,2.0);
00069                         
00070                 }
00071                 densities.at(i) = densities.at(i-1)* Random.GetNumber(1.1,1.5);
00072                 resistivities.at(i) = Random.GetNumber(0.0,3.0);
00073                 velocities.at(i) = velocities.at(i-1)* Random.GetNumber(1.05,1.2);
00074         }
00075         transform(resistivities.begin(),resistivities.end(),resistivities.begin(),boost::bind(fopow<double,double>(),10,_1));
00076         MTSynth.SetResistivities(resistivities);
00077         MTSynth.SetThicknesses(mtthickness);
00078         MTSynth.writemodel(outputbase+"_mt.mod");
00079         MTSynth.writeplot(outputbase+"_mt.plot");
00080         MTSynth.GetData();
00081         MTSynth.WriteAsMtt(outputbase);
00082         
00083         RecModel.Init(Configuration.nlayers); //initialize model
00084         RecModel.slowness = Configuration.slowness; //copy slowness
00085         copy(recthickness.begin(),recthickness.end(),RecModel.Thickness.begin()); //copy values
00086         copy(velocities.begin(),velocities.end(),RecModel.SVelocity.begin());
00087         copy(densities.begin(),densities.end(),RecModel.Density.begin());
00088         transform(RecModel.SVelocity.begin(),RecModel.SVelocity.end(), //PVelocity is calculated by multiplication of SVel by poisson's ratio
00089                         RecModel.PVelocity.begin(),boost::bind(multiplies<double>(),_1,Configuration.poisson));
00090         RecModel.dt = Configuration.dt; // set dt
00091         RecModel.npts = Configuration.npts; //set number of points
00092         RecSynth.CalcRecSynth(Configuration.outputbase+"_rec",RecModel,RecSynthData,false); // Calculate forward model
00093         RecModel.WritePlot(outputbase+"_rec");
00094         RecSynthData.WriteAsSac(outputbase+".rec");
00095 }

Generated on Fri Jul 4 15:30:20 2008 for GPLIB++ by  doxygen 1.5.5