enumerate.cpp

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <algorithm>
#include <numeric>
#include <boost/bind.hpp>
#include <boost/shared_ptr.hpp>
#include "Iso1DMTObjective.h"
#include "MTStation.h"
#include "C1DRecObjective.h"
#include "SeismicDataComp.h"
#include "gentypes.h"
#include "convert.h"
#include "NetCDFTools.h"
#include "VecMat.h"
#include "SeisTools.h"
using namespace std;




int main()
{
        ofstream outfile;
        MTStation MTData;


        const int nlayers = 2;
        ttranscribed mtvalues(nlayers *2 ),recvalues(nlayers *2);
        double mtresult, recresult;

        MTData.GetData("nobias.mtt");
        boost::shared_ptr<const SeismicDataComp> RecData(new SeismicDataComp("smp.rec", SeismicDataComp::sac));
         Normalize(const_cast<SeismicDataComp*>(RecData.get())->GetData());

        outfile.open("misfit.2d");
        Iso1DMTObjective MTObjective(MTData);
        MTObjective.AppendFitParameters(&MTTensor::GetRhoxy, &MTTensor::GetdRhoxy, 0.05);

        C1DRecObjective RecObjective(RecData,0, 4.0, 0.001, 0.07);
        RecObjective.SetPoisson(sqrt(3.0));
        RecObjective.SetTimeWindow(1.0,30.0);
        RecObjective.SetErrorLevel(0.02);

        mtvalues(0) = 2;
        mtvalues(1) = 1.3;
        mtvalues(2) = 20;
        mtvalues(3) = 10;
        recvalues(0) = 20;
        recvalues(1) = 0;
        recvalues(2) = 3.46;
        recvalues(3) = 4.6;
        double jbase = 25;
        double velbase = 3.0;
        double resbase = 1.5;
        string filename;
        const int thicksteps = 50;
        const int velsteps = 100;

        gplib::rmat mtfit(thicksteps,velsteps), recfit(thicksteps,velsteps);
        gplib::rvec thicknesses(thicksteps), velocities(velsteps), resistivities(velsteps);
        ofstream recthick("rec.thick");
        for (int j = 0; j < thicksteps; ++j)
        {
                        double currthick =jbase + 0.5*j;

                        mtvalues(2) = currthick;
                        recvalues(0) = currthick;
                        thicknesses(j) =currthick;
                        cout << "J: " << j << endl;
                        for (int k = 0; k < velsteps; ++k)
                        {
                                double currres = resbase + 0.01*k;
                                double currvel = velbase + 0.01*k;
                                velocities(k) = currvel;
                                resistivities(k) = currres;
                                mtvalues(0) = currres;
                                recvalues(2) = currvel;
                                //cout << currvel << endl;
                                filename = "out"+stringify(currthick)+stringify(currvel);
                                recresult = RecObjective.CalcPerformance(recvalues);
                                mtresult = MTObjective.CalcPerformance(mtvalues);
                                recfit(j,k) = recresult;
                                mtfit(j,k) = mtresult;
                                RecObjective.WriteData(filename+".rec");
                                RecObjective.WriteModel(filename+".mod");
                        }
                        outfile << endl << endl;
        }
        Write2DDataAsNetCDF("mt.nc",mtfit, thicknesses, resistivities,"Thickness", "Resistivity", "Misfit");
        Write2DDataAsNetCDF("rec.nc",recfit, thicknesses, velocities,"Thickness", "Velocity", "Misfit");
        outfile.close();
}