intlf.cpp

Go to the documentation of this file.

#include <iostream>
#include <string>
#include "ThreeDMTModel.h"
#include <fstream>

using namespace std;
/*! /file
 * Plot the integrated conductivity of a 3D Model up to a certain depth.
 */
int main(void)
{
        C3DModel *Model;
        string infilename,outfilename;
        double maxdepth;
        const double kmlatitude = 111194;
        const double kmlongitude = 71474;       
        /*const double xmin = 48;
        const double ymin = 2.;*/
        //const double xmin = 44.5;
        //const double ymin = -1.;
        double xmin, ymin;
        int mode;
                
        cout << "Modelfile: ";
        cin >> infilename;
        cout << "Type 1 for old Mackie, 2 for Winglink: ";
        cin  >> mode;
        cout << "Depth: ";
        cin >> maxdepth;
        cout << "Xmin in degree: ";
        cin >> xmin;
        cout << "Ymin in degree: ";
        cin >> ymin;
        cout << "Outfile: ";
        cin >> outfilename;
        if (mode == 1)
                Model = new 3DMackieModel;
        else
                Model = new C3DWinglink;
        Model->GetData(infilename);
        Array<double,2> SumCond(Model->xsize.size(),Model->ysize.size());
        Array<double,1> xscale(Model->xsize.size()), yscale(Model->ysize.size());
        SumCond = 0;
        double currentdepth = 0;
        int layerindex = 0;
        while (currentdepth < maxdepth)
        {       
                for (int i = 0; i < Model->xsize.size(); ++i)
                        for (int j = 0; j < Model->ysize.size(); ++j)
                                SumCond(i,j) += 1./Model->Resistivities(Model->Model(i,j,layerindex)) * Model->zsize(layerindex);
                currentdepth += Model->zsize(layerindex)/1000;
                cout << currentdepth << endl;
                ++layerindex;
        }
        xscale = 0;
        yscale = 0;
        for (int i = 1; i < Model->xsize.size(); ++i)
                xscale(i) = xscale(i-1) + Model->xsize(i-1);
        for (int i = 1; i < Model->ysize.size(); ++i)
                yscale(i) = yscale(i-1) + Model->ysize(i-1);
        ofstream outfile(outfilename.c_str());
        for (int i = 0; i < Model->xsize.size(); ++i)
                        for (int j = 0; j < Model->ysize.size(); ++j)
                        {
                                outfile << ymin + yscale(j)/kmlongitude << " " << xmin + xscale(i)/kmlatitude << " " << SumCond(i,j) << endl; 
                        }
}