C1DAnisoMTSynthData.cpp

Go to the documentation of this file.

#include "C1DAnisoMTSynthData.h"
#include <cassert>
#include <limits>
#include <fstream>
#include "FatalException.h"

namespace gplib
  {
    const int frequenzen = 50;
    const float T[frequenzen] =
      { 0.0025, 0.00308337, 0.00513901, 0.00625, 0.0104167, 0.0125, 0.0208333,
          0.025, 0.0416667, 0.05, 0.0833333, 0.1, 0.166667, 0.2, 0.333333, 0.4,
          0.667, 0.64, 1.06667, 1.29333, 2.15554, 2.59997, 4.33332, 5.19994,
          8.66701, 10.666667, 16.000, 21.333, 32.000, 42.667, 64.0, 85.332,
          127.992, 170.678, 256.016, 341.297, 512.033, 682.594, 1023.541,
          1366.120, 2049.180, 2732.240, 4098.361, 5464.481, 8196.722,
          10869.565, 16393.443, 28800, 43200, 86400 };

    void C1DAnisoMTSynthData::TransformRho()
      {
        assert(rho1.size() == rho2.size());
        assert(rho1.size() == rho3.size());
        assert(rho1.size() == strikes.size());
        assert(rho1.size() == slants.size());
        assert(rho1.size() == dips.size());

        const unsigned int nlayers = rho1.size();
        effcond1.assign(nlayers, 0.0);
        effcond2.assign(nlayers, 0.0);
        effstrike.assign(nlayers, 0.0);
        for (unsigned int i = 0; i < nlayers; ++i)
          {
            double s1 = 1. / rho1.at(i);//transform to conductivity
            double s2 = 1. / rho2.at(i);
            double s3 = 1. / rho3.at(i);
            double strike = PI * strikes.at(i) / 180.0; //transform angles to radian
            double dip = PI * dips.at(i) / 180.0;
            double slant = PI * slants.at(i) / 180.0;
            //calculate some frequently used sines and cosines
            double sps = sin(strike);
            double cps = cos(strike);
            double sth = sin(dip);
            double cth = cos(dip);
            double sfi = sin(slant);
            double cfi = cos(slant);
            double pom1 = s1 * cfi * cfi + s2 * sfi * sfi;
            double pom2 = s1 * sfi * sfi + s2 * cfi * cfi;
            double pom3 = (s1 - s2) * sfi * cfi;
            double c2ps = cps * cps;
            double s2ps = sps * sps;
            double c2th = cth * cth;
            double s2th = sth * sth;
            double csps = cps * sps;
            double csth = cth * sth;
            //contruct the conductivity tensor elements
            double sg11 = pom1 * c2ps + pom2 * s2ps * c2th - 2. * pom3 * cth
                * csps + s3 * s2th * s2ps;
            double sg12 = pom1 * csps - pom2 * c2th * csps + pom3 * cth * (c2ps
                - s2ps) - s3 * s2th * csps;
            double sg13 = -pom2 * csth * sps + pom3 * sth * cps + s3 * csth
                * sps;
            double sg21 = sg12;
            double sg22 = pom1 * s2ps + pom2 * c2ps * c2th + 2. * pom3 * cth
                * csps + s3 * s2th * c2ps;
            double sg23 = pom2 * csth * cps + pom3 * sth * sps - s3 * csth
                * cps;
            double sg31 = sg13;
            double sg32 = sg23;
            double sg33 = pom2 * s2th + s3 * c2th;
            double axx = sg11 - sg13 * sg31 / sg33;
            double axy = sg12 - sg13 * sg32 / sg33;
            double ayx = sg21 - sg31 * sg23 / sg33;
            double ayy = sg22 - sg23 * sg32 / sg33;
            //calculate effective horizontal conductivity
            double da12 = sqrt((axx - ayy) * (axx - ayy) + 4.0 * axy * ayx);
            effcond1.at(i) = (0.5 * (axx + ayy + da12));
            effcond2.at(i) = (0.5 * (axx + ayy - da12));
            if (gsl_fcmp(da12, 0.0, GSL_DBL_EPSILON) != 0)
              {
                effstrike.at(i) = 0.5 * acos((axx - ayy) / da12);
                assert(gsl_finite(effstrike.at(i)));
              }
            else
              {
                effstrike.at(i) = 0.0;
              }
            if (axy < 0.0)
              effstrike.at(i) *= -1.0;
            //std::cout << effcond1.at(i) << " " << effcond2.at(i) << " " << effstrike.at(i)*180./PI << std::endl;

          }
      }

    void C1DAnisoMTSynthData::GetData()
      {

        calc_frequencies = GetFrequencies();
        if (calc_frequencies.empty())
          {
            for (int i = 0; i < frequenzen; ++i)
              calc_frequencies.push_back(1.0 / T[i]);
          }
        Assign(calc_frequencies.size());
        TransformRho();
        CalcZ();
        Update();

      }

    inline dcomp dfp(const dcomp x)
      {
        //std::cout << 1.0 + exp(-2.0 * x) << std::endl;
        return 1.0 + exp(-2.0 * x);
      }

    inline dcomp dfm(const dcomp x)
      {
        //std::cout << 1.0 - exp(-2.0 * x) << std::endl;
        return 1.0 - exp(-2.0 * x);
      }

    void C1DAnisoMTSynthData::CalcZ()
      {
        dcomp omegamu;
        const double convfactor = 1. / (1000. * mu);

        for (unsigned int i = 0; i < calc_frequencies.size(); ++i)
          {
            omegamu = -I * 8e-7 * PI * PI * calc_frequencies.at(i);
            dcomp k0 = (1. - I) * 2.0e-3 * PI * sqrt(calc_frequencies.at(i)
                / 10.0);
            dcomp zxyrot = k0 / sqrt(effcond1.back());
            dcomp zyxrot = -k0 / sqrt(effcond2.back());
            double currstrike = effstrike.back();
            MTTensor CurrZ(dcomp(0.0, 0.0), zxyrot, zyxrot, dcomp(0.0, 0.0),
                calc_frequencies.at(i));
            for (int j = thicknesses.size() - 2; j >= 0; --j)
              {
                double currthick = 1000.0 * thicknesses.at(j);
                if (currstrike != effstrike.at(j) && effcond1.at(j)
                    != effcond2.at(j))
                  {
                    CurrZ.Rotate(effstrike.at(j) - currstrike);
                    currstrike = effstrike.at(j);
                  }
                MTTensor BottomZ(CurrZ);
                //std::cout << CurrZ.GetZxy() * convfactor<< " " << CurrZ.GetZyx() * convfactor << std::endl;
                dcomp dz1 = k0 / sqrt(effcond1.at(j));
                dcomp dz2 = k0 / sqrt(effcond2.at(j));
                dcomp ag1 = k0 * sqrt(effcond1.at(j)) * currthick;
                dcomp ag2 = k0 * sqrt(effcond2.at(j)) * currthick;

                dcomp detzbottom = BottomZ.GetDet();
                dcomp denominator = detzbottom * dfm(ag1) * dfm(ag2) / (dz1
                    * dz2) + BottomZ.GetZxy() * dfm(ag1) * dfp(ag2) / dz1
                    - BottomZ.GetZyx() * dfp(ag1) * dfm(ag2) / dz2 + dfp(ag1)
                    * dfp(ag2);
                //std::cout << dz1 << dz2 << ag1 << ag2 << detzbottom << denominator << std::endl;
                CurrZ.SetZxx() = 4.0 * BottomZ.GetZxx() * exp(-ag1 - ag2)
                    / denominator;
                CurrZ.SetZxy() = (BottomZ.GetZxy() * dfp(ag1) * dfp(ag2)
                    - BottomZ.GetZyx() * dfm(ag1) * dfm(ag2) * dz1 / dz2
                    + detzbottom * dfp(ag1) * dfm(ag2) / dz2 + dfm(ag1) * dfp(
                    ag2) * dz1) / denominator;
                CurrZ.SetZyx() = (BottomZ.GetZyx() * dfp(ag1) * dfp(ag2)
                    - BottomZ.GetZxy() * dfm(ag1) * dfm(ag2) * dz2 / dz1
                    - detzbottom * dfm(ag1) * dfp(ag2) / dz1 - dfp(ag1) * dfm(
                    ag2) * dz2) / denominator;
                CurrZ.SetZyy() = 4.0 * BottomZ.GetZyy() * exp(-ag1 - ag2)
                    / denominator;
              }
            CurrZ.SetRotangle() = currstrike;
            CurrZ.SetZxx() = convfactor * conj(CurrZ.SetZxx());
            CurrZ.SetZxy() = convfactor * conj(CurrZ.SetZxy());
            CurrZ.SetZyx() = convfactor * conj(CurrZ.SetZyx());
            CurrZ.SetZyy() = convfactor * conj(CurrZ.SetZyy());
            assert(gsl_finite(abs(CurrZ.GetZxx()))); //there seem to be some problems with certain calculations
            assert(gsl_finite(abs(CurrZ.GetZxy()))); //safeguard to help investigate
            assert(gsl_finite(abs(CurrZ.GetZyx())));
            assert(gsl_finite(abs(CurrZ.GetZyy())));
            //std::cout << CurrZ.GetZxy() << " " << CurrZ.GetZyx() << std::endl << std::endl;
            CurrZ.Rotate(-CurrZ.GetRotangle());
            //
            SetMTData().at(i) = CurrZ;
            /*std::copy(effstrike.begin(),effstrike.end(),std::ostream_iterator<double>(std::cout," "));
             std::cout << std::endl;
             std::copy(effcond1.begin(),effcond1.end(),std::ostream_iterator<double>(std::cout," "));
             std::cout << std::endl;
             std::copy(effcond2.begin(),effcond2.end(),std::ostream_iterator<double>(std::cout," "));
             std::cout << std::endl;*/

          }
      }

    void C1DAnisoMTSynthData::WriteModel(std::string filename)
      {
        assert(rho1.size() == thicknesses.size());
        assert(rho2.size() == thicknesses.size());
        assert(rho3.size() == thicknesses.size());
        assert(strikes.size() == thicknesses.size());
        assert(slants.size() == thicknesses.size());
        assert(dips.size() == thicknesses.size());
        std::ofstream outfile(filename.c_str());
        if (outfile)
          {
            outfile << thicknesses.size() << std::endl;
            for (unsigned int i = 0; i < thicknesses.size(); ++i)
              {
                outfile << thicknesses.at(i) << "  " << rho1.at(i) << " ";
                outfile << rho2.at(i) << "  " << rho3.at(i) << " ";
                outfile << strikes.at(i) << "  " << slants.at(i) << " ";
                outfile << dips.at(i) << std::endl;
              }
          }
        else
          {
            throw FatalException("Problem writing to file: " + filename);
          }

      }

    void C1DAnisoMTSynthData::ReadModel(std::string filename)
      {

        int layers;
        double currthick, currho1, currho2, currho3, currstrike, currslant,
            currdip;
        int i = 0;

        std::ifstream infile(filename.c_str());
        infile >> layers;
        if (!infile.good() || layers < 0)
          throw FatalException("Problem reading file: " + filename);
        thicknesses.reserve(layers);
        rho1.reserve(layers);
        rho2.reserve(layers);
        rho3.reserve(layers);
        strikes.reserve(layers);
        slants.reserve(layers);
        dips.reserve(layers);
        while (infile.good())
          {
            infile >> currthick >> currho1 >> currho2 >> currho3 >> currstrike
                >> currslant >> currdip;
            if (infile.good())
              {
                thicknesses.push_back(currthick);
                rho1.push_back(currho1);
                rho2.push_back(currho2);
                rho3.push_back(currho3);
                strikes.push_back(currstrike);
                slants.push_back(currslant);
                dips.push_back(currdip);
                ++i;
              }
          }
        assert(rho1.size() == thicknesses.size());
        assert(rho2.size() == thicknesses.size());
        assert(rho3.size() == thicknesses.size());
        assert(strikes.size() == thicknesses.size());
        assert(slants.size() == thicknesses.size());
        assert(dips.size() == thicknesses.size());
        assert(i==layers);
      }

    void C1DAnisoMTSynthData::WritePlot(std::string filename)
      {
        std::ofstream outfile(filename.c_str());
        double cumthick = 0;
        TransformRho();
        assert(effcond1.size() == thicknesses.size());
        assert(effcond2.size() == thicknesses.size());
        assert(effstrike.size() == thicknesses.size());
        for (unsigned int i = 0; i < thicknesses.size(); ++i)
          {
            outfile << cumthick << "  " << 1. / effcond1.at(i) << " ";
            outfile << 1. / effcond2.at(i) << "  " << effstrike.at(i) * 180.
                / PI;
            outfile << std::endl;
            cumthick += thicknesses.at(i);
            outfile << cumthick << "  " << 1. / effcond1.at(i) << " ";
            outfile << 1. / effcond2.at(i) << "  " << effstrike.at(i) * 180.
                / PI;
            outfile << std::endl;
          }
      }

    gplib::rvec C1DAnisoMTSynthData::GetModelVector()
      {
        TransformRho();
        assert(effcond1.size() == thicknesses.size());
        assert(effcond2.size() == thicknesses.size());
        assert(effstrike.size() == thicknesses.size());
        const unsigned int length = thicknesses.size();
        gplib::rvec result(length * 4);
        for (unsigned int i = 0; i < length; ++i)
          {
            result(i) = log10(1. / effcond1.at(i));
            result(i + length) = thicknesses.at(i);
            result(i + 2 * length) = log10(effcond1.at(i) / effcond2.at(i));
            result(i + 3 * length) = effstrike.at(i) * 180. / PI;
          }
        return result;
      }

    C1DAnisoMTSynthData::C1DAnisoMTSynthData(const C1DAnisoMTSynthData &old) :
      MTStation(old), calc_frequencies(old.calc_frequencies), thicknesses(
          old.thicknesses), strikes(old.strikes), slants(old.slants), dips(
          old.dips), rho1(old.rho1), rho2(old.rho2), rho3(old.rho3), effcond1(
          old.effcond1), effcond2(old.effcond2), effstrike(old.effstrike)
      {

      }
    C1DAnisoMTSynthData::C1DAnisoMTSynthData()
      {
      }

    C1DAnisoMTSynthData::~C1DAnisoMTSynthData()
      {
      }
  }