testrecobjective.cpp

Go to the documentation of this file.

#include <iostream>
#include "AbsVelRecObjective.h"
#include "SurfaceWaveData.h"
#include "ResPkModel.h"
#include "SeismicDataComp.h"
#include "CCalcRecConf.h"
#include <string>
#include <boost/bind.hpp>
#include "SeisTools.h"
using namespace std;
int main()
{
        string datafilename, modelfilename;
        cout << "Base for Datafiles: ";
        cin >> datafilename;
        cout << "Name of Modelfile: ";
        cin >> modelfilename;
        double slowness;
        cout << "Slowness: ";
        cin >> slowness;
        CCalcRecConf Config;
        Config.GetData("calcrec.conf");
        SurfaceWaveData SurfData;
        boost::shared_ptr<const SeismicDataComp> RecData(new SeismicDataComp(datafilename+".rec", SeismicDataComp::sac));
         Normalize(const_cast<SeismicDataComp*>(RecData.get())->GetData());
        SurfData.ReadAscii(datafilename+".vel");
        ResPkModel Model;
        Model.GetData(modelfilename);

        AbsVelRecObjective SeisObjective(RecData,SurfData,Config.shift,Config.sigma,Config.cc,slowness,RecCalc::iterdecon);

        SeisObjective.SetFitExponent(2);
        SeisObjective.SetPoisson(1.7320508);
        SeisObjective.SetErrorLevel(0.02);
        SeisObjective.SetAbsVelWeight(1);
        SeisObjective.SetRecWeight(0);
        SeisObjective.SetTimeWindow(0,40);
        const unsigned int nlayers = Model.GetThickness().size();
        const unsigned int nparams = nlayers * 2;
        ttranscribed member(nparams);
        for (unsigned i = 0; i < nlayers; ++i)
        {
                member(i) = Model.SetThickness().at(i);
                member(i+nlayers) = Model.SetSVelocity().at(i);
        }
        cout << "Misfit: " << SeisObjective.CalcPerformance(member) << endl;
        SeisObjective.WriteData("obj");
        SeisObjective.WriteModel("obj.mod");
}

---