All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines
Public Member Functions
QGpCoreWave::MagnetoTelluricFactory Class Reference

Brief description of class still missing. More...

#include <MagnetoTelluricFactory.h>

List of all members.

Public Member Functions

const QVector< double > & angularFrequencies () const
bool calculate (const Resistivity1DModel &model, double staticShift=1.0)
MagnetoTelluricCurve curve () const
const ComplexValueimpedanceXY () const
void linkX (QList< MagnetoTelluricCurve > &curves) const
void linkX (MagnetoTelluricCurve &c) const
 MagnetoTelluricFactory ()
void readReport (QDataStream &s)
void setAngularFrequency ()
void setX (const QList< MagnetoTelluricCurve > &curves)
void setX (const MagnetoTelluricCurve &c)
void setX (const QVector< double > &x)
void toStream (QTextStream &sOut, const MagnetoTelluricPointOptions &options) const
void writeReport (QDataStream &s) const
 ~MagnetoTelluricFactory ()

Detailed Description

Brief description of class still missing.

Full description of class still missing


Constructor & Destructor Documentation

References TRACE.

  {
    TRACE;
    _impedanceXY=0;
  }

References TRACE.

  {
    TRACE;
    delete [] _impedanceXY;
  }

Member Function Documentation

bool QGpCoreWave::MagnetoTelluricFactory::calculate ( const Resistivity1DModel model,
double  staticShift = 1.0 
)

Cagniard algorithm written by Max Moorkamp, ported to this library by Marc Wathelet.

Permission to include this code here must be requested from Max (GPL into LGPL library-add defines). Easier to include that code here than to link to gplib which makes reference to many other libraries.

References QGpCoreTools::conjugate(), QGpCoreWave::Resistivity1DModel::depths(), QGpCoreTools::exp(), QGpCoreWave::Resistivity1DModel::layerCount(), MAGNETIC_CONSTANT, QGpCoreTools::Complex::re, QGpCoreWave::Resistivity1DModel::resistivities(), QGpCoreTools::Complex::setIm(), QGpCoreTools::Value< numberType >::setValue(), and QGpCoreTools::sqrt().

Referenced by DinverDCCore::TargetList::magnetoTelluricMisfit(), main(), and MagnetoTelluricReader::parse().

  {
    if(model.layerCount()==0) {
      return false;
    }
    Complex omegamu, omega;
    Complex kcurr, klow;
    Complex xi;
    QVector<Complex> alpha;
    Complex adm;
    double d;
    double sigmacurr, sigmalow;
    int nx=_x.count();
    int nLayers=model.layerCount();
    const QVector<double>& depths=model.depths();
    const QVector<double>& resistivities=model.resistivities();
    for(int i=0; i<nx; i++) {
      omega.setIm(-_x.at(i));
      omegamu.setIm(-MAGNETIC_CONSTANT*_x.at(i));
      //printf("omegamu[%i]=%lf+i*%lf\n", i, omegamu.re(), omegamu.im());
      alpha.fill(Complex::null, resistivities.count());
      if(nLayers<2) {
        sigmacurr=1.0/resistivities.last();
        kcurr=sqrt(omegamu * sigmacurr);
        //printf("kcurr(HS)=%lf+i*%lf\n",kcurr.re(), kcurr.im());
      } else {
        for(int j=nLayers-2; j>=0; --j) {
          //printf("Layer %i, low %i\n", j, j+1);
          sigmalow=1.0/resistivities.at(j + 1);
          sigmacurr=1.0/resistivities.at(j);
          kcurr=sqrt(omegamu * sigmacurr);
          klow=sqrt(omegamu * sigmalow);
          //printf("kcurr=%lf+i*%lf\n",kcurr.re(), kcurr.im());
          //printf("klow=%lf+i*%lf\n",klow.re(), klow.im());
          if(kcurr.re() < 0.0) {
            kcurr *= -1.0;
            klow *= -1.0;
          }
          Complex km=kcurr + klow;
          km*=km;
          xi=omegamu * (sigmacurr - sigmalow)/km;
          //printf("xi=%lg+i*%lg\n",xi.re(),xi.im());
          alpha[j + 1]=(xi + alpha.at(j + 1))/(1.0+xi * alpha.at(j + 1));
          //printf("alpha[j+1]=%lg+i*%lg\n",alpha.at(j+1).re(),alpha.at(j+1).im());
          d=j>0 ? depths.at(j)-depths.at(j-1) : depths.at(j);
          //printf("d=%lf\n",d);
          alpha[j]=alpha.at(j + 1) * exp(-2. * kcurr * d);
          //printf("alpha[j]=%lg+i*%lg\n",alpha.at(j).re(),alpha.at(j).im());
        }
      }
      //printf("kcurr(surf)=%lf+i*%lf\n",kcurr.re(), kcurr.im());
      adm=kcurr/omega * ((1.0- alpha.at(0))/(1.0 + alpha.at(0)));
      //printf("adm=%lg+i*%lg\n",adm.re(), adm.im());
      _impedanceXY[i].setValue(staticShift*conjugate(1./adm));
      //printf("Zxy=%lg+i*%lg\n",_impedanceXY[i].value().re(), _impedanceXY[i].value().im());
    }
    return true;
  }

References QGpCoreTools::Curve< pointType >::resize(), QGpCoreTools::StatisticalValue< numberType >::setMean(), QGpCoreTools::StatisticalPoint< numberType >::setX(), and TRACE.

  {
    TRACE;
    MagnetoTelluricCurve c;
    int n=_x.count();
    c.resize(n);
    double factor=0.5/M_PI;
    for(int i=0; i<n; i++) {
      MagnetoTelluricPoint& p=c[i];
      p.setX(factor*_x[i]);
      p.setMean(_impedanceXY[i].value());
    }
    return c;
  }

Repeat calls to linkX()

References TRACE.

Referenced by main(), and DinverDCCore::TargetList::validateTargets().

  {
    TRACE;
    for(QList<MagnetoTelluricCurve>::iterator it=curves.begin();it!=curves.end();++it) {
      linkX(*it);
    }
  }

References QGpCoreTools::endl(), QGpCoreWave::MagnetoTelluricCurve::linkX(), QGpCoreTools::tr(), and TRACE.

  {
    TRACE;
    if(_x.isEmpty()) {
      App::stream() << tr("### ERROR ### : frequency list not set or empty when linking frequencies.") << endl;
      return;
    }
    c.linkX(_x);
  }

References QGpCoreTools::Value< numberType >::setValue(), and TRACE.

Referenced by DinverDCGui::MagnetoTelluricViewer::report2plot(), and DinverDCCore::DCReportBlock::write().

  {
    TRACE;
    int nFrequencies;
    s >> nFrequencies;
    _x.resize(nFrequencies);
    for(int i=0; i<nFrequencies; i++) {
      s >> _x[i];
    }
    Complex c;
    delete _impedanceXY;
    _impedanceXY=new ComplexValue[nFrequencies];
    for(int i=0; i<nFrequencies; i++) {
      s >> c;
      _impedanceXY[i].setValue(c);
    }
  }

Switch to omega from frequency. Must be called only once.

References TRACE.

Referenced by main(), MagnetoTelluricReader::parse(), and DinverDCCore::TargetList::validateTargets().

  {
    TRACE;
    double factor=2*M_PI;
    for(QVector<double>::iterator it=_x.begin(); it!=_x.end(); it++ ) {
      (*it) *= factor;;
    }
  }

Repeat calls to setX()

References TRACE.

Referenced by main(), MagnetoTelluricReader::parse(), and DinverDCCore::TargetList::validateTargets().

  {
    TRACE;
    for(QList<MagnetoTelluricCurve>::const_iterator it=curves.begin();it!=curves.end();++it) {
      setX(*it);
    }
  }

References QGpCoreTools::Curve< pointType >::begin(), QGpCoreTools::Curve< pointType >::count(), QGpCoreTools::Curve< pointType >::end(), TRACE, and QGpCoreTools::unique().

  {
    TRACE;
    _x.reserve(_x.count()+c.count());
    for(MagnetoTelluricCurve::const_iterator it=c.begin(); it!=c.end(); it++ ) {
      _x.append(it->x());
    }
    qSort(_x);
    unique(_x);
    // Allocates values for storing calculated impedances
    delete [] _impedanceXY;
    _impedanceXY=new ComplexValue[_x.count()];
  }
void QGpCoreWave::MagnetoTelluricFactory::setX ( const QVector< double > &  x)

References TRACE, and QGpCoreTools::unique().

  {
    TRACE;
    _x << x;
    qSort(_x);
    unique(_x);
    // Allocates values for storing calculated impedances
    delete [] _impedanceXY;
    _impedanceXY=new ComplexValue[_x.count()];
  }
void QGpCoreWave::MagnetoTelluricFactory::toStream ( QTextStream &  sOut,
const MagnetoTelluricPointOptions options 
) const

References QGpCoreTools::flush(), QGpCoreTools::Point2D::setX(), QGpCoreTools::Point2D::setY(), QGpCoreWave::MagnetoTelluricPointOptions::toDouble(), TRACE, QGpCoreTools::Value< numberType >::value(), and QGpCoreTools::Point2D::x().

Referenced by main(), and MagnetoTelluricReader::parse().

  {
    TRACE;
    double factor=0.5/M_PI;
    int n=_x.count();
    Point2D p;
    for(int i=0; i<n; i++) {
      p.setX(factor*_x[i]);
      p.setY(options.toDouble(p.x(), _impedanceXY[i].value()));
      sOut << p << "\n";
    }
    sOut << flush;
  }
void QGpCoreWave::MagnetoTelluricFactory::writeReport ( QDataStream &  s) const

References TRACE, and QGpCoreTools::Value< numberType >::value().

Referenced by DinverDCCore::DCReportBlock::write().

  {
    TRACE;
    int nFrequencies=_x.count();
    s << nFrequencies;
    for(int i=0; i<nFrequencies; i++) {
      s << _x.at(i);
    }
    for(int i=0; i<nFrequencies; i++) {
      s << _impedanceXY[i].value();
    }
  }

The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines