MTStation.cpp

Go to the documentation of this file.

#include "MTStation.h"
#include "miscfunc.h"
//we only include the parsers for .j and .edi files
//if we have a proper installation of antlr
#ifdef HAVEANTLR
#include "EDILexer.hpp"
#include "EDIParser.hpp"
#include "JLexer.hpp"
#include "JParser.hpp"
#endif
#include <iomanip>
#include <iterator>
#include <string>
#include "FatalException.h"
#include <boost/filesystem/operations.hpp>
#include <boost/filesystem/convenience.hpp>
#include <boost/bind.hpp>
#include <functional>
#include <fstream>
#include <cassert>
#include "Util.h"
using namespace std;

namespace gplib
  {
    void TrimFilename(string &name)
      {
        size_t slashpos = name.find('/', 0);
        if (slashpos != string::npos)
          name = name.substr(slashpos + 1);
        size_t dotpos = name.find('.', 0);
        if (dotpos != string::npos)
          name.erase(dotpos);
      }

    void MTStation::InitialSetup()
      {
        latitude = 0;
        longitude = 0;
        elevation = 0;
        azimuth = 0;
      }

    MTStation::MTStation()
      {
        InitialSetup();
      }

    MTStation::~MTStation()
      {
      }

    MTStation::MTStation(const std::string filename)
      {
        InitialSetup();
        GetData(filename);
      }

    MTStation::MTStation(const int size) :
      MTData(size)
      {
        InitialSetup();
        Assign(size);
      }
    void MTStation::AssignAll(const int nfreq)
      {
        Assign(nfreq);
      }
    void MTStation::Assign(const int nfreq)
      {
        MTData.assign(nfreq, MTTensor());
        TFData.assign(nfreq, MagneticTF());
      }

    /*! The Update() method first updates all components individually and then
     * recalculates all derived quantities that depend on more than one component
     */

    void MTStation::Update()
      {

      }

    /*! Rotates the component indicated by i by angle rotangle. Does not Update. */
    void MTStation::DoRotate(const int i, const double angle)
      {
        dcomp newxx, newxy, newyx, newyy;
        newxx = MTData.at(i).Zxx * std::pow(cos(angle), 2) + (MTData.at(i).Zxy
            + MTData.at(i).Zyx) * sin(angle) * cos(angle) + MTData.at(i).Zyy
            * std::pow(sin(angle), 2);
        newxy = MTData.at(i).Zxy * std::pow(cos(angle), 2) - (MTData.at(i).Zxx
            - MTData.at(i).Zyy) * sin(angle) * cos(angle) - MTData.at(i).Zyx
            * std::pow(sin(angle), 2);
        newyx = MTData.at(i).Zyx * std::pow(cos(angle), 2) - (MTData.at(i).Zxx
            - MTData.at(i).Zyy) * sin(angle) * cos(angle) - MTData.at(i).Zxy
            * std::pow(sin(angle), 2);
        newyy = MTData.at(i).Zyy * std::pow(cos(angle), 2) - (MTData.at(i).Zxy
            + MTData.at(i).Zyx) * sin(angle) * cos(angle) + MTData.at(i).Zxx
            * std::pow(sin(angle), 2);
        MTData.at(i).Zxx = newxx;
        MTData.at(i).Zxy = newxy;
        MTData.at(i).Zyx = newyx;
        MTData.at(i).Zyy = newyy;
        MTData.at(i).rotangle += angle;
      }

    /*! Rotate(const double rotangle) rotates the impedance tensor by the angle rotangle in radian
     * and updates all derived quantities
     */
    void MTStation::Rotate(const double rotangle)
      {
        for (unsigned int i = 0; i < MTData.size(); ++i)
          {
            DoRotate(i, rotangle);
          }
        Update();
      }

    /*! Rotate(void) rotates the impedance tensor by the negative of the angle given in the rotangles field
     */
    void MTStation::Rotate(void)
      {
        for (unsigned int i = 0; i < MTData.size(); ++i)
          {
            DoRotate(i, -MTData.at(i).rotangle);
          }
        Update();
      }
    /*! Return a vector containing the available frequencies; This should be rewritten to be more efficient
     * */
    trealdata MTStation::GetFrequencies() const
      {
        trealdata temp(MTData.size());
        transform(MTData.begin(), MTData.end(), temp.begin(), mem_fun_ref(
            &MTTensor::GetFrequency));
        return temp;
      }

    void MTStation::SetFrequencies(const trealdata &freqs)
      {
        //assert(MTData.size() == freqs.size());
        const unsigned int nfreqs = freqs.size();
        if (MTData.size() != nfreqs) //this invalidates the data stored there before
          {
            MTData.resize(nfreqs);
            TFData.resize(nfreqs);
          }
        for (unsigned int i = 0; i < nfreqs; ++i)
          MTData.at(i).frequency = freqs.at(i);
        assert(MTData.size() == TFData.size());

      }
#ifdef HAVEANTLR
    void MTStation::ReadEdi(const std::string filename)
      {
        using namespace antlr;
        ifstream infile(filename.c_str());
        const double invalid = 1e30;

        if (infile)
          {
            try
              {
                EDILexer lexer(infile);
                EDIParser parser(lexer);
                parser.edi_file();
                latitude = parser.latitude;
                longitude = parser.longitude;
                elevation = parser.elevation;
                azimuth = parser.azimuth;
                name = filename;
                TrimFilename(name);
                unsigned int valfreqs = 0;
                for (unsigned int i = 0; i < parser.frequency.size(); ++i)
                  {
                    if (abs(parser.DataXX.at(i)) > invalid)
                      {
                        parser.frequency.at(i) = 0;
                      }
                    else
                      {
                        valfreqs++;
                      }
                  }
                //frequency.assign(parser.frequency.size(),0);
                //copy(parser.frequency.begin(),parser.frequency.end(),frequency.begin());
                //rotangles.assign(parser.frequency.size(),0);
                azimuth += parser.rotangles.at(0) * 180. / PI;
                Assign(valfreqs);
                unsigned int mtindex = 0;
                //There might be some invalid entries that we do not want to copy, we therefore
                //go through all parsed entries with the for loop, but have a separate index
                //for our storage class
                for (unsigned int i = 0; i < parser.rotangles.size(); ++i)
                  {
                    if (abs(parser.DataXX.at(i)) < invalid)
                      {
                        MTData.at(mtindex).frequency = parser.frequency.at(i);
                        MTData.at(mtindex).rotangle = parser.rotangles.at(i)
                        * PI / 180.;
                        MTData.at(mtindex).Zxx = parser.DataXX.at(i);
                        MTData.at(mtindex).Zxy = parser.DataXY.at(i);
                        MTData.at(mtindex).Zyx = parser.DataYX.at(i);
                        MTData.at(mtindex).Zyy = parser.DataYY.at(i);
                        MTData.at(mtindex).dZxx = parser.dDataXX.at(i);
                        MTData.at(mtindex).dZxy = parser.dDataXY.at(i);
                        MTData.at(mtindex).dZyx = parser.dDataYX.at(i);
                        MTData.at(mtindex).dZyy = parser.dDataYY.at(i);
                        mtindex++;
                      }
                  }
                mtindex = 0;
                if (parser.DataZY.empty() == false)
                  {
                    for (unsigned int i = 0; i < parser.DataZY.size(); ++i)
                      {
                        if (abs(parser.DataZX.at(i)) < invalid)
                          {
                            TFData.at(mtindex).frequency = parser.frequency.at(
                                i);
                            TFData.at(mtindex).Tx = parser.DataZX.at(i);
                            TFData.at(mtindex).Ty = parser.DataZY.at(i);
                            TFData.at(mtindex).dTx = parser.dDataZX.at(i);
                            TFData.at(mtindex).dTy = parser.dDataZY.at(i);
                            mtindex++;
                          }
                      }
                  }
              }
            catch (ANTLRException& e)
              {
                cerr << "Parse exception: " << e.toString() << endl;
              }
            Update();
          }
        else
          {
            throw FatalException("File not found: " + filename);
          }
      }
#endif

    void MTStation::ReadPek1D(const std::string filename)
      {
        ifstream infile(filename.c_str());
        if (infile)
          {
            const double convfactor = 1. / (1000. * mu);
            double period, rxx, imxx, rxy, imxy, ryx, imyx, ryy, imyy;
            while (infile.good())
              {
                infile >> period >> rxx >> imxx >> rxy >> imxy >> ryx >> imyx
                    >> ryy >> imyy;

                if (infile.good())
                  {
                    MTTensor CurrentImp;
                    CurrentImp.frequency = 1. / period;
                    CurrentImp.Zxx = (rxx + I * imxx) * convfactor;
                    CurrentImp.Zxy = (rxy + I * imxy) * convfactor;
                    CurrentImp.Zyx = (ryx + I * imyx) * convfactor;
                    CurrentImp.Zyy = (ryy + I * imyy) * convfactor;
                    MTData.push_back(CurrentImp);
                    TFData.push_back(MagneticTF());
                  }
              }
            name = filename;
            TrimFilename(name);
          }
        else
          {
            throw FatalException("File not found: " + filename);
          }
      }

#ifdef HAVEANTLR
    void MTStation::ReadJ(const std::string filename)
      {

        using namespace antlr;
        ifstream infile(filename.c_str());

        if (infile)
          {
            try
              {
                JLexer lexer(infile);
                JParser parser(lexer);
                parser.jfile();
                latitude = parser.latitude;
                longitude = parser.longitude;
                elevation = parser.elevation;
                azimuth = parser.azimuth;
                name = parser.name;
                TrimFilename(name);
                Assign(parser.frequency.size());

                //frequency = parser.frequency;
                if (parser.zassigned == false && parser.rassigned == false)
                  {
                    cerr << "No MT data in file !" << endl;
                  }
                else
                  {
                    for (unsigned int i = 0; i < parser.frequency.size(); ++i)
                      {
                        MTData.at(i).frequency = parser.frequency.at(i);
                        MTData.at(i).rotangle = parser.azimuth;
                        MTData.at(i).Zxx = parser.DataXX.at(i);
                        MTData.at(i).Zxy = parser.DataXY.at(i);
                        MTData.at(i).Zyx = parser.DataYX.at(i);
                        MTData.at(i).Zyy = parser.DataYY.at(i);
                        MTData.at(i).dZxx = parser.dDataXX.at(i);
                        MTData.at(i).dZxy = parser.dDataXY.at(i);
                        MTData.at(i).dZyx = parser.dDataYX.at(i);
                        MTData.at(i).dZyy = parser.dDataYY.at(i);
                        MTData.at(i).Rx = parser.Rx.at(i);
                        MTData.at(i).Ry = parser.Ry.at(i);
                      }
                  }
                if (parser.tassigned)
                  {
                    for (unsigned int i = 0; i < parser.DataZY.size(); ++i)
                      {
                        TFData.at(i).frequency = parser.frequency.at(i);
                        TFData.at(i).Tx = parser.DataZX.at(i);
                        TFData.at(i).Ty = parser.DataZY.at(i);
                        TFData.at(i).dTx = parser.dDataZX.at(i);
                        TFData.at(i).dTy = parser.dDataZY.at(i);
                        TFData.at(i).Rz = parser.Rz.at(i);
                      }
                  }
              }
            catch (ANTLRException& e)
              {
                cerr << "Parse exception: " << e.toString() << endl;
              }
            Update();
          }
        else
          {
            throw FatalException("File not found: " + filename);
          }

      }
#endif

    void MTStation::ReadMtt(const std::string filename)
      {
        ifstream infile;
        double currentreal, currentimag;
        int nentries = 0;
        infile.open(filename.c_str());
        while (infile.good())
          {
            infile >> currentreal;
            ++nentries;
          }
        infile.close();
        if (((nentries - 1) % 23) != 0)
          throw FatalException("Number of records does not match expected: "
              + filename);
        const int nrecords = (nentries - 1) / 23;
        Assign(nrecords);
        infile.open(filename.c_str());
        int currentrecord = 0;
        if (infile)
          {
            while (infile.good()) // read in inputfile
              {
                infile >> currentreal;
                if (infile.good())
                  {
                    MTData.at(currentrecord).frequency = currentreal;
                    TFData.at(currentrecord).frequency = currentreal;
                    infile >> currentreal;
                    MTData.at(currentrecord).Nu = currentreal;
                    infile >> currentreal >> currentimag;
                    MTData.at(currentrecord).Zxx = (currentreal + I
                        * currentimag);
                    infile >> currentreal >> currentimag;
                    MTData.at(currentrecord).Zxy = (currentreal + I
                        * currentimag);
                    infile >> currentreal >> currentimag;
                    MTData.at(currentrecord).Zyx = (currentreal + I
                        * currentimag);
                    infile >> currentreal >> currentimag;
                    MTData.at(currentrecord).Zyy = (currentreal + I
                        * currentimag);
                    infile >> currentreal;
                    MTData.at(currentrecord).dZxx = currentreal;
                    infile >> currentreal;
                    MTData.at(currentrecord).dZxy = currentreal;
                    infile >> currentreal;
                    MTData.at(currentrecord).dZyx = currentreal;
                    infile >> currentreal;
                    MTData.at(currentrecord).dZyy = currentreal;
                    infile >> currentreal >> currentimag;
                    TFData.at(currentrecord).Tx = currentreal + I * currentimag;
                    infile >> currentreal >> currentimag;
                    TFData.at(currentrecord).Ty = currentreal + I * currentimag;
                    infile >> currentreal;
                    TFData.at(currentrecord).dTx = currentreal;
                    infile >> currentreal;
                    TFData.at(currentrecord).dTy = currentreal;
                    infile >> currentreal;
                    MTData.at(currentrecord).Rx = currentreal;
                    infile >> currentreal;
                    MTData.at(currentrecord).Ry = currentreal;
                    infile >> currentreal;
                    TFData.at(currentrecord).Rz = currentreal;
                    ++currentrecord;
                  }
              }
            infile.close();
            name = filename;
            TrimFilename(name);
          }
        else
          {
            throw FatalException("File not found: " + filename);
          }
        Update();
      }

    void MTStation::GetData(const string filename)
      {
        if (!boost::filesystem::exists(filename))
          throw FatalException("File does not exist : " + filename);
        string ending = boost::filesystem::extension(filename);
        if (ending == "")
          {
            cerr << "File has no extension ! " << filename << endl;
            dataformat = unknown;
            throw FatalException("File has no extension ! " + filename);
          }
        if (ending == ".mtt")
          {
            ReadMtt(filename);
            dataformat = mtt;
          }
        else if (ending == ".dat" || ending == ".j")
          {
#ifdef HAVEANTLR
            ReadJ(filename);
            dataformat = j;
#else
            throw FatalException("Did not compile with .j support");
#endif
          }
        else if (ending == ".edi")
          {
#ifdef HAVEANTLR
            ReadEdi(filename);
            dataformat = edi;
#else
            throw FatalException("Did not compile with .edi support");
#endif
          }
        else if (ending == ".pek")
          {
            ReadPek1D(filename);
            dataformat = pek;
          }
        else
          {
            dataformat = unknown;
            throw FatalException("File not recognised: " + filename);
          }
        assert(MTData.size() == TFData.size());
        Rotate(); //rotate data to 0 degree

      }

    void MTStation::WriteData(const std::string filename)
      {
        name = filename;
        WriteBack();
      }

    void MTStation::WriteBack()
      {
        switch (dataformat)
          {
        case mtt:
          WriteAsMtt(name.c_str());
          break;
        case j:
          WriteAsJ(name.c_str());
          break;
        case edi:
          WriteAsEdi(name.c_str());
          break;
        default:
          throw FatalException(
              "Unknown File Format for Writing ! This should not happen");
          }

      }

    void MTStation::WriteAsMtt(const string filename)
      {
        Update();
        WriteMtt((filename + ".mtt").c_str());
      }

    void MTStation::WriteAsEdi(const string filename)
      {
        throw FatalException("WriteAsEdi not implemented yet !");
      }

    void MTStation::WriteJBlock(boost::function<complex<double> (
        const MTTensor*)> Comp, boost::function<double(const MTTensor*)> Err,
        ofstream &outfile, const double convfactor)
      {
        for (unsigned int i = 0; i < MTData.size(); ++i)
          {
            if (MTData.at(i).frequency != 0)
              {
                outfile << setfill(' ') << setw(15)
                    << resetiosflags(ios::fixed) << 1. / MTData.at(i).frequency
                    << " ";
                outfile << setfill(' ') << setw(15) << convfactor
                    * boost::bind(Comp, &MTData.at(i))().real() << " ";
                outfile << setfill(' ') << setw(15) << convfactor
                    * boost::bind(Comp, &MTData.at(i))().imag() << " ";
                outfile << setfill(' ') << setw(15) << convfactor
                    * boost::bind(Err, &MTData.at(i))() << " ";
                outfile << setfill(' ') << setw(15) << setiosflags(ios::fixed)
                    << 1.000 << " " << endl;
              }
          }
      }

    void MTStation::WriteAsJ(const string filename)
      {
        ofstream outfile;
        outfile.open((filename + ".j").c_str());
        int actfreqs = 0;
        const double convfactor = 4. * acos(-1.) * 1e-4;

        outfile << "# J-File Produced by CJData.cpp" << endl;
        outfile << ">LATITUDE  = " << latitude << endl;
        outfile << ">LONGITUDE = " << longitude << endl;
        outfile << ">ELEVATION = " << elevation << endl;
        outfile << ">AZIMUTH   = " << azimuth << endl;
        outfile << name << endl;
        for (unsigned int i = 0; i < MTData.size(); ++i)
          if (MTData.at(i).frequency != 0)
            actfreqs++;
        outfile << "ZXX S.I." << endl;
        outfile << actfreqs << endl;
        //WriteJBlock(DataXX,outfile,convfactor);
        WriteJBlock(&MTTensor::GetZxx, &MTTensor::GetdZxx, outfile, convfactor);
        outfile << "ZXY S.I." << endl;
        outfile << actfreqs << endl;
        WriteJBlock(&MTTensor::GetZxy, &MTTensor::GetdZxy, outfile, convfactor);

        outfile << "ZYX S.I." << endl;
        outfile << actfreqs << endl;
        WriteJBlock(&MTTensor::GetZyx, &MTTensor::GetdZyx, outfile, convfactor);

        outfile << "ZYY S.I." << endl;
        outfile << actfreqs << endl;
        WriteJBlock(&MTTensor::GetZyy, &MTTensor::GetdZyy, outfile, convfactor);

        /*      if (norm(DataZY.Z.at(0)) !=0)
         {
         outfile << "TZX " << endl;
         outfile << actfreqs << endl;
         WriteJBlock(DataZX,outfile,1.);
         }

         if (norm(DataZY.Z.at(0)) !=0)
         {
         outfile << "TZY " << endl;
         outfile << actfreqs << endl;
         WriteJBlock(DataZY,outfile,1.);
         }*/
        outfile.close();
      }

    void MTStation::WriteMtt(const std::string filename)
      {
        ofstream outfile;

        outfile.open(filename.c_str());
        vector<int> mttindex(MTData.size());
        // index(frequency, mttindex);
        assert(MTData.size() == TFData.size());
        for (unsigned int i = 0; i < MTData.size(); ++i) //write mtt-file
          {
            //k = mttindex.at(i);
            if (MTData.at(i).frequency > 0)
              {
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::scientific) << MTData.at(i).frequency;
                outfile << "  " << resetiosflags(ios::scientific)
                    << setprecision(5) << MTData.at(i).Nu << "\n";

                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).Zxx.real()
                    << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).Zxx.imag()
                    << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).Zxy.real()
                    << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).Zxy.imag()
                    << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).Zyx.real()
                    << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).Zyx.imag()
                    << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).Zyy.real()
                    << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).Zyy.imag()
                    << " ";
                outfile << "\n";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).dZxx << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).dZxy << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).dZyx << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).dZyy << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << TFData.at(i).Tx.real() << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << TFData.at(i).Tx.imag() << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << TFData.at(i).Ty.real() << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << TFData.at(i).Ty.imag() << " ";
                outfile << "\n";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << TFData.at(i).dTx << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << TFData.at(i).dTy << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).Rx << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << MTData.at(i).Ry << " ";
                outfile << setw(9) << setfill(' ') << setprecision(4)
                    << setiosflags(ios::fixed) << TFData.at(i).Rz << " ";
                outfile << "\n" << resetiosflags(ios::fixed);
                //write out mtt file entries
              }

          }
        outfile.close();
      }
    MTStation::MTStation(const MTStation &old) :
      latitude(old.latitude), longitude(old.longitude),
          elevation(old.elevation), name(old.name), azimuth(old.azimuth),
          MTData(old.MTData), TFData(old.TFData),

          dataformat(old.dataformat)
      {
        Update();
      }
    MTStation& MTStation::operator=(const MTStation& source)
      {
        if (this == &source)
          return *this;

        copy(source.MTData.begin(), source.MTData.end(), back_inserter(
            this->MTData));
        copy(source.TFData.begin(), source.TFData.end(), back_inserter(
            this->TFData));

        this->latitude = source.latitude;
        this->longitude = source.longitude;
        this->elevation = source.elevation;
        this->azimuth = source.azimuth;
        this->name = source.name;
        this->dataformat = source.dataformat;
        this->Update();
        return *this;
      }
  }