makeinput.cpp

Go to the documentation of this file.

#include <iostream>
#include "C1DMTSynthData.h"
#include "CResPkModel.h"
#include "CRecCalc.h"
#include "CUniformRNG.h"
#include "MakeInputConf.h"
#include "CSeismicDataComp.h"
#include "Adaptors.h"
#include <vector>
#include <boost/bind.hpp>
#include <string>

using namespace std;
string version = "$Id: makeinput.cpp 1093 2007-05-08 17:28:32Z max $";

int main(int argc, char *argv[])
{
        cout << " This is makeinput: Create MT and receiver function models and data" << endl; // write some info
        cout << " to be used as input for inversion methods." << endl;
        cout << " If one argument is given, this argument is taken as the filename base for all output files." << endl;
        cout << " Otherwise the output file base is read from the configuration file 'makeinput.conf' " << endl;
        cout << " The configuration file contains all other relevant settings. " << endl;
        cout << " This is Version: " << version << endl << endl;
        UniformRNG Random;
        MakeInputConf Configuration;
        Configuration.GetData("makeinput.conf");
        
        string outputbase;
        if (argc > 1)
                outputbase = argv[1];
        else
                outputbase = Configuration.outputbase;
                
        C1DMTSynthData MTSynth;
        ResPkModel RecModel;
        RecCalc RecSynth(Configuration.shift, Configuration.omega, Configuration.sigma, Configuration.cc);
        SeismicDataComp RecSynthData;
        const double basevel = 3.2;
        const double basedens = 2.6;
        const double baseres = 1.0 + Random.GetNumber();
        const double basethick = 1.0 + Random.GetNumber()*10;
        vector<double> velocities(Configuration.nlayers);
        vector<double> densities(Configuration.nlayers);
        vector<double> resistivities(Configuration.nlayers);
        vector<double> mtthickness(Configuration.nlayers);
        vector<double> recthickness(Configuration.nlayers);
        densities.front() = basedens;
        resistivities.front() = baseres;
        velocities.front() = basevel;
        mtthickness.front() = basethick;
        if (Configuration.correlated)
        {
                recthickness.front() = basethick;
        }
        else
        {
                recthickness.front() = basethick * Random.GetNumber(1.1,2.0);
        }
        for (int i = 1; i < Configuration.nlayers; ++i)
        {
                mtthickness.at(i) = mtthickness.at(i-1) * Random.GetNumber(1.1,2.0);
                if (Configuration.correlated)
                {
                         recthickness.at(i) = mtthickness.at(i);
                }
                else
                {
                                recthickness.at(i) = mtthickness.at(i) * Random.GetNumber(1.1,2.0);
                        
                }
                densities.at(i) = densities.at(i-1)* Random.GetNumber(1.1,1.5);
                resistivities.at(i) = Random.GetNumber(0.0,3.0);
                velocities.at(i) = velocities.at(i-1)* Random.GetNumber(1.05,1.2);
        }
        transform(resistivities.begin(),resistivities.end(),resistivities.begin(),boost::bind(fopow<double,double>(),10,_1));
        MTSynth.SetResistivities(resistivities);
        MTSynth.SetThicknesses(mtthickness);
        MTSynth.writemodel(outputbase+"_mt.mod");
        MTSynth.writeplot(outputbase+"_mt.plot");
        MTSynth.GetData();
        MTSynth.WriteAsMtt(outputbase);
        
        RecModel.Init(Configuration.nlayers); //initialize model
        RecModel.slowness = Configuration.slowness; //copy slowness
        copy(recthickness.begin(),recthickness.end(),RecModel.Thickness.begin()); //copy values
        copy(velocities.begin(),velocities.end(),RecModel.SVelocity.begin());
        copy(densities.begin(),densities.end(),RecModel.Density.begin());
        transform(RecModel.SVelocity.begin(),RecModel.SVelocity.end(), //PVelocity is calculated by multiplication of SVel by poisson's ratio
                        RecModel.PVelocity.begin(),boost::bind(multiplies<double>(),_1,Configuration.poisson));
        RecModel.dt = Configuration.dt; // set dt
        RecModel.npts = Configuration.npts; //set number of points
        RecSynth.CalcRecSynth(Configuration.outputbase+"_rec",RecModel,RecSynthData,false); // Calculate forward model
        RecModel.WritePlot(outputbase+"_rec");
        RecSynthData.WriteAsSac(outputbase+".rec");
}

---