A layered model with tilted interfaces. More...

#include <RefractionDippingModel.h>

Inheritance diagram for QGpCoreWave::RefractionDippingModel:

Classes
class	RefractionContext
class	RefractionStorage
class	VariableDL
class	VariableDR
class	VariableN
class	VariableV
Public Types
enum	Direction { LeftToRight, RightToLeft }
Public Member Functions
void	begin ()
virtual GeophysicalModel *	clone () const
double	depthLeft (int iLayer) const
double	depthRight (int iLayer) const
void	end ()
virtual GeophysicalContext *	expressionContext () const
void	fromSeismic1DModel (const Seismic1DModel &ref, const QVector< double > &pitch, Seismic1DModel::SlownessType slowType)
virtual bool	fromStream (QTextStream &s, QString *comments=0)
void	fromStream (QDataStream &s)
bool	isCriticalLayer (Direction dir, int iLayer) const
virtual bool	isEmpty () const
int	layerCount () const
double	propagate (Direction dir, int iLayer, int iDeepestLayer, double &x) const
Curve< Point2D >	ray (const TiltNode &src, const TiltNode &rec)
	RefractionDippingModel ()
	RefractionDippingModel (int nLayers)
	RefractionDippingModel (const RefractionDippingModel &o)
double	setDepthLeft (int iLayer, double d)
double	setDepthRight (int iLayer, double d)
double	setSlowness (int iLayer, double s)
void	setXLeft (double x)
void	setXRight (double x)
double	slowness (int iLayer) const
virtual void	toStream (QTextStream &s) const
void	toStream (QDataStream &s) const
virtual QString	toString () const
double	travelTime (const TiltNode &src, const TiltNode &rec, int &minTimeLayer)
double	xLeft () const
double	xRight () const
	~RefractionDippingModel ()
Static Public Member Functions
static QString	formatHelp ()

Detailed Description

A layered model with tilted interfaces.

Each layer has one unique slowness. The bottom interface geometry is defined by its left and right depths. Depths are defined at xLeft() and xRight().

Member Enumeration Documentation

enum QGpCoreWave::RefractionDippingModel::Direction

Enumerator:

LeftToRight
RightToLeft

{LeftToRight, RightToLeft};

Constructor & Destructor Documentation

QGpCoreWave::RefractionDippingModel::RefractionDippingModel ( )

References TRACE.

{
  TRACE;
  initMembers();
}

QGpCoreWave::RefractionDippingModel::RefractionDippingModel ( int nLayers )

References TRACE.

{
  TRACE;
  _layerCount=nLayers;
  initMembers();
  allocateData();
}

QGpCoreWave::RefractionDippingModel::RefractionDippingModel ( const RefractionDippingModel & o )

References TRACE.

  : GeophysicalModel()
{
  TRACE;
  _layerCount=o._layerCount;
  initMembers();
  allocateData();
  for(int i=1;i<_layerCount;i++) {
    _depthL[i]=o._depthL[i];
    _depthR[i]=o._depthR[i];
    _slow[i]=o._slow[i];
    _minX[i]=o._minX[i];
    _maxX[i]=o._maxX[i];
  }
}

QGpCoreWave::RefractionDippingModel::~RefractionDippingModel ( )

References end(), and TRACE.

{
  TRACE;
  if(_cosangle) end();
  deleteData();
}

Member Function Documentation

void QGpCoreWave::RefractionDippingModel::begin ( )

Geometrical and physical parameters are supposed to be all initialized.

Compute relative angles between layers. For each layer, the critical angle angle is computed. Incidence angle are then back propagated to the surface (refraction at angle lower than critical).

References QGpCoreTools::asin(), QGpCoreTools::cos(), QGpCoreTools::endl(), QGpCoreTools::sin(), QGpCoreTools::tr(), and TRACE.

Referenced by QGpCompatibility::CompatHodoChroneData::calculate(), QGpCoreWave::RefractionFactory::calculate(), and RefraReader::parse().

{
  TRACE;
  _cosangle=new double [_layerCount];
  _sinangle=new double [_layerCount];
  _critl=new double *[_layerCount];
  _critr=new double *[_layerCount];
  _shift=new double [_layerCount];
  _h=new double [_layerCount];
  double * angle=new double [_layerCount];

  // Init computed parameter for first layer
  angle[0]=atan((_depthL[0]-_depthR[0])/(_xR-_xL));
  _cosangle[0]=cos(angle[0]);
  _sinangle[0]=sin(angle[0]);
  _minX[0]=0.5*(_xR+_xL);
  _maxX[0]=0.5*(_xR+_xL);
  _critl[0]=0;
  _critr[0]=0;
  for(int i=1;i<_layerCount;i++) {
    angle[i]=atan((_depthL[i]-_depthR[i])/(_xR-_xL));
    _cosangle[i]=cos(angle[i]);
    _sinangle[i]=sin(angle[i]);
    double tmp=(_depthL[i]-_depthL[i-1])/cos(angle[i]-angle[i-1]);
    _shift[i]=-tmp*_sinangle[i-1];
    _h[i]=tmp*_cosangle[i];
    _minX[i]=0.5*(_xR+_xL);
    _maxX[i]=0.5*(_xR+_xL);
    if(_slow[i]<_slow[i-1]) { // refraction at critical angle
      double inc0=asin(_slow[i]/_slow[i-1]);
      _critl[i]=new double [i];
      _critl[i][i-1]=inc0;
      int j;
      for(j=i-2;j>=0;j--) {
        double inc=_critl[i][j+1]-(angle[j+2]-angle[j+1]);
        double a=_slow[j+1]/_slow[j]*sin(inc);
        if(a<1.0) _critl[i][j]=asin(a);
        else {
          App::stream() << tr("Critical angle cannot be reached in layer %1 from LeftToRight").arg(j) << endl;
          delete [] _critl[i]; // this is not a good candidate for critical angle
          _critl[i]=0;
        }
      }
      _critr[i]=new double [i];
      _critr[i][i-1]=-inc0;
      for(j=i-2;j>=0;j--) {
        double inc=_critr[i][j+1]-(angle[j+2]-angle[j+1]);
        double a=_slow[j+1]/_slow[j]*sin(inc);
        if(a<1.0) _critr[i][j]=asin(a);
        else {
          App::stream() << tr("Critical angle cannot be reached in layer %1 from RightToLeft").arg(j) << endl;
          delete [] _critr[i]; // this is not a good candidate for critical angle
          _critr[i]=0;
        }
      }
    } else {
      _critl[i]=0;
      _critr[i]=0;
    }
  }
  delete [] angle;
}

virtual GeophysicalModel* QGpCoreWave::RefractionDippingModel::clone ( ) const [inline, virtual]

Implements QGpCoreWave::GeophysicalModel.

{return new RefractionDippingModel(*this);}

double QGpCoreWave::RefractionDippingModel::depthLeft ( int iLayer ) const [inline]

Referenced by QGpCoreWave::RefractionDippingModel::VariableDL::value().

{return _depthL[iLayer];}

double QGpCoreWave::RefractionDippingModel::depthRight ( int iLayer ) const [inline]

Referenced by QGpCoreWave::RefractionDippingModel::VariableDR::value().

{return _depthR[iLayer];}

void QGpCoreWave::RefractionDippingModel::end ( )

References TRACE.

Referenced by QGpCompatibility::CompatHodoChroneData::calculate(), QGpCoreWave::RefractionFactory::calculate(), RefraReader::parse(), and ~RefractionDippingModel().

{
  TRACE;
  delete [] _cosangle;
  delete [] _sinangle;
  delete [] _shift;
  delete [] _h;
  _cosangle=0;
  _sinangle=0;
  _shift=0;
  _h=0;
  int i;
  for(i=1;i<_layerCount;i++) delete [] _critl[i];
  delete [] _critl;
  for(i=1;i<_layerCount;i++) delete [] _critr[i];
  delete [] _critr;
  _critl=0;
  _critr=0;
}

virtual GeophysicalContext* QGpCoreWave::RefractionDippingModel::expressionContext ( ) const [inline, virtual]

Implements QGpCoreWave::GeophysicalModel.

{return new RefractionContext;}

QString QGpCoreWave::RefractionDippingModel::formatHelp ( ) [static]

References TRACE.

Referenced by QGpCoreWave::GeophysicalModel::allFormatHelp().

{
  TRACE;
  return "  Line 1    <number of layers including half-space for first model>\n"
         "  Line 2    <x left depth profile> <x right depth profile>\n"
         "  Line 3    <left thickness (m)> <right thickness (m)> <V (m/s)>\n"
         "  ....\n"
         "  Line n    0 0 <V (m/s)>\n"
         "  Line n+1  <number of layers including half-space for second model>\n"
         "  ....";
}

void QGpCoreWave::RefractionDippingModel::fromSeismic1DModel	(	const Seismic1DModel &	ref,
		const QVector< double > &	pitch,
		Seismic1DModel::SlownessType	slowType
	)

Convert ref horizontal layered model to a titl layered model using pitch profile. The layer pitch is a real number between -1 and 1. pitch profile must contain exactly the same number of layers than ref.

-1 left thickness is forced to 0
0 means no pitch.
1 right thickness is forced to 0

References QGpCoreWave::Seismic1DModel::h(), QGpCoreWave::Seismic1DModel::layerCount(), QGpCoreWave::Seismic1DModel::P, QGpCoreWave::Seismic1DModel::S, QGpCoreWave::Seismic1DModel::slowP(), QGpCoreWave::Seismic1DModel::slowS(), and TRACE.

Referenced by DinverDCCore::TargetList::refractionMisfit().

{
  TRACE;
  if(_layerCount!=ref.layerCount()) {
    _layerCount=ref.layerCount();
    deleteData();
    allocateData();
  }
  ASSERT(_layerCount==pitch.count());
  int n1=_layerCount-1;
  for(int i=0; i<n1;i++) {
    switch(slowType) {
    case Seismic1DModel::P:
      _slow[i]=ref.slowP(i);
      break;
    case Seismic1DModel::S:
      _slow[i]=ref.slowS(i);
      break;
    }
    _depthL[i+1]=_depthL[i] + (1.0+pitch[i])*ref.h(i);
    _depthR[i+1]=_depthR[i] + (1.0-pitch[i])*ref.h(i);
  }
  switch(slowType) {
  case Seismic1DModel::P:
    _slow[n1]=ref.slowP(n1);
    break;
  case Seismic1DModel::S:
    _slow[n1]=ref.slowS(n1);
    break;
  }
}

bool QGpCoreWave::RefractionDippingModel::fromStream	(	QTextStream &	s,
		QString *	comments = `0`
	)		`[virtual]`

Loads model from a text stream

Format:

      layerCount
      xLeft xRight
      thicknessLeft thicknessRight V (layer 1)
      ...
      0 0 V (layer n, half-space)

Returns false if an error occured during reading or if the read model is not consistent. It returns true even if the number of layers is null or cannot be read. This is the responsability of the caller to test if the number of layers is greater than 1 before doing anything else.

Implements QGpCoreWave::GeophysicalModel.

References QGpCoreTools::endl(), QGpCoreTools::StringSection::isEmpty(), QGpCoreTools::StringSection::isValid(), QGpCoreTools::StringSection::nextField(), QGpCoreTools::StringSection::set(), QGpCoreTools::StringSection::toDouble(), QGpCoreTools::StringSection::toInt(), QGpCoreTools::StringSection::toString(), QGpCoreTools::tr(), TRACE, and QGpCoreTools::StringSection::trimmed().

Referenced by RefraReader::parse(), and refractionModelMode().

{
  TRACE;
  // read the number of layers
  StringSection field, fields;
  const QChar * ptr;
  QString buf;
  buf=File::readLineNoComments(s, comments);
  fields.set(buf);
  fields.trimmed();
  ptr=0;
  field=fields.nextField(ptr);
  if(field.isValid()) {
    bool ok=true;
    _layerCount=field.toInt(&ok);
    if(!ok && !field.isEmpty()) { // Not a blank line but failed
      App::stream() << tr("Tilted layered model: bad number of layers: \"%1\" interpreted as %2").arg(field.toString()).arg(_layerCount) << endl;
       _layerCount=0;
       return false;
    }
  } else {
    _layerCount=0;
    return true;  // This is just a blank line, before doing anything, you must test if the number of layer is >0
  }
  if(_layerCount>0 && !s.atEnd()) {
    allocateData();
    buf=File::readLineNoComments(s, comments);
    fields.set(buf);
    fields.trimmed();
    ptr=0;
    field=fields.nextField(ptr);
    if(field.isValid()) {
      _xL=field.toDouble();
    } else return false;
    field=fields.nextField(ptr);
    if(field.isValid()) {
      _xR=field.toDouble();
    } else return false;
    // read the data
    _depthL[0]=0.0;
    _depthR[0]=0.0;
    for(int i=0;i<_layerCount;i++) {
      if(s.atEnd()) {
        App::stream() << tr("Tilted layered model: reaching the end of file, %1 layers found, expecting %2")
                  .arg(i).arg(_layerCount) << endl;
        _layerCount=0;
        return false;
      }
      buf=File::readLineNoComments(s, comments);
      fields.set(buf);
      fields.trimmed();
      ptr=0;
      field=fields.nextField(ptr);
      if(i<_layerCount-1) {
        if(!setValue(_depthL[i+1], field, i, tr("left thickness"), false) ) return false;
        _depthL[i+1] += _depthL[i];
      }
      field=fields.nextField(ptr);
      if(i<_layerCount-1) {
        if(!setValue(_depthR[i+1], field, i, tr("right thickness"), false) ) return false;
        _depthR[i+1] += _depthR[i];
     }
      field=fields.nextField(ptr);
      if(!setValue(_slow[i], field, i, tr("V"), false) ) return false;
      // Switch to slowness
      _slow[i]=1/_slow[i];
    }
    return true;
  } else {
    _layerCount=0;
    return false;
  }
}

void QGpCoreWave::RefractionDippingModel::fromStream ( QDataStream & s )

bool QGpCoreWave::RefractionDippingModel::isCriticalLayer	(	Direction	dir,
		int	iLayer
	)		const `[inline]`

References LeftToRight, and TRACE.

{
  TRACE;
  switch(dir) {
  case LeftToRight:
    return _critl[iLayer];
    break;
  default:
    return _critr[iLayer];
    break;
  }
}

virtual bool QGpCoreWave::RefractionDippingModel::isEmpty ( ) const [inline, virtual]

Implements QGpCoreWave::GeophysicalModel.

{return layerCount()==0;}

int QGpCoreWave::RefractionDippingModel::layerCount ( ) const [inline]

Referenced by RefraReader::parse(), refractionModelMode(), QGpCoreWave::RefractionDippingModel::VariableDL::setValue(), QGpCoreWave::RefractionDippingModel::VariableDR::setValue(), QGpCoreWave::RefractionDippingModel::VariableV::setValue(), QGpCoreWave::RefractionDippingModel::VariableDL::value(), QGpCoreWave::RefractionDippingModel::VariableDR::value(), and QGpCoreWave::RefractionDippingModel::VariableV::value().

{return _layerCount;}

double QGpCoreWave::RefractionDippingModel::propagate	(	Direction	dir,
		int	iLayer,
		int	iDeepestLayer,
		double &	x
	)		const

Return time needed to propagate across layer iLayer if the target layer is iDeepestLayer. x is set the local coordinate of the ray at the bottom interface. Local coordinates are counted from the reference profile.

References QGpCoreTools::cos(), LeftToRight, RightToLeft, SAFE_UNINITIALIZED, QGpCoreTools::sin(), QGpCoreTools::sqrt(), and TRACE.

{
  TRACE;
  double inc;
  SAFE_UNINITIALIZED(inc,0);
  switch(dir) {
  case LeftToRight:
    inc=_critl[iDeepestLayer][iLayer-1];
    break;
  case RightToLeft:
    inc=_critr[iDeepestLayer][iLayer-1];
    break;
  }
  double cosab, cosb, sinb, d, dt;
  d=sqrt(x*x+_h[iLayer]*_h[iLayer]);
  cosb=_h[iLayer]/d;
  sinb=x/d;
  cosab=_cosangle[iLayer]*cosb-_sinangle[iLayer]*sinb; // cos(a+b)
  double cosinc=cos(inc);
  double dcosinc=d/cosinc;
  // time to cross the layer
  dt=_slow[iLayer-1]*cosab*dcosinc;
  /* intersection of the ray with interface counted from the reference profile
     sin(a+b-i)=sin(a+b)cos(i)-cos(a+b)sin(i)
               =[sin(a)cos(b)+cos(a)sin(b)]cos(i)-cos(a+b)sin(i) */
  x=((_sinangle[iLayer]*cosb+_cosangle[iLayer]*sinb)*cosinc-cosab*sin(inc))*
    dcosinc-_shift[iLayer];
  return dt;
}

Curve< Point2D > QGpCoreWave::RefractionDippingModel::ray	(	const TiltNode &	src,
		const TiltNode &	rec
	)

Return ray path between two nodes src and rec.

References QGpCoreWave::TiltPath::abscissa(), QGpCoreTools::Curve< pointType >::resize(), and travelTime().

Referenced by RefraReader::parse().

{
  Curve<Point2D> rayPoints;
  const TiltPath * srcPaths, * recPaths;
  if(!node2paths(src, rec, srcPaths, recPaths)) return rayPoints;
  int endLayer;
  travelTime(srcPaths, recPaths, endLayer);
  double x;
  int i;
  int lastpray=endLayer*2+1;
  rayPoints.resize(lastpray+1);
  int pathIndex=endLayer-1;
  if(pathIndex<0) pathIndex=0;

  for(i=0;i<=endLayer;i++) {
    x=srcPaths[pathIndex].abscissa(i);
    rayPoints[i].setX(_xL+x*_cosangle[i] );
    rayPoints[i].setY(_depthL[i]-x*_sinangle[i] );
  }

  for(i=0;i<=endLayer;i++) {
    x=recPaths[pathIndex].abscissa(i);
    rayPoints[lastpray-i].setX(_xL+x*_cosangle[i] );
    rayPoints[lastpray-i].setY(_depthL[i]-x*_sinangle[i] );
  }

  return rayPoints;
}

double QGpCoreWave::RefractionDippingModel::setDepthLeft	(	int	iLayer,
		double	d
	)		`[inline]`

Referenced by QGpCoreWave::RefractionDippingModel::VariableDL::setValue().

{return _depthL[iLayer]=d;}

double QGpCoreWave::RefractionDippingModel::setDepthRight	(	int	iLayer,
		double	d
	)		`[inline]`

Referenced by QGpCoreWave::RefractionDippingModel::VariableDR::setValue().

{return _depthR[iLayer]=d;}

double QGpCoreWave::RefractionDippingModel::setSlowness	(	int	iLayer,
		double	s
	)		`[inline]`

Referenced by QGpCoreWave::RefractionDippingModel::VariableV::setValue().

{return _slow[iLayer]=s;}

void QGpCoreWave::RefractionDippingModel::setXLeft ( double x ) [inline]

Referenced by QGpCompatibility::CompatHodoChroneData::calculate(), QGpCoreWave::RefractionFactory::calculate(), and outputDCModel().

{_xL=x;}

void QGpCoreWave::RefractionDippingModel::setXRight ( double x ) [inline]

Referenced by QGpCompatibility::CompatHodoChroneData::calculate(), and QGpCoreWave::RefractionFactory::calculate().

{_xR=x;}

double QGpCoreWave::RefractionDippingModel::slowness ( int iLayer ) const [inline]

Referenced by QGpCoreWave::RefractionDippingModel::VariableV::value().

{return _slow[iLayer];}

void QGpCoreWave::RefractionDippingModel::toStream ( QTextStream & s ) const [virtual]

Saves model to a text stream

Format:

      layerCount
      xLeft xRight
      thicknessLeft thicknessRight V (layer 1)
      ...
      0 0 V (layer n, half-space)

Implements QGpCoreWave::GeophysicalModel.

References QGpCoreTools::flush(), and TRACE.

Referenced by refractionModelMode().

{
  TRACE;
  static const QString str3("%1 %2 %3\n");
  s << _layerCount << "\n"
    << _xL << " " << _xR << "\n";
  for(int i=0;i<_layerCount-1;i++)
    s << str3.arg(_depthL[i+1]-_depthL[i], 0, 'g', 20)
             .arg(_depthR[i+1]-_depthR[i], 0, 'g', 20)
             .arg(1/_slow[i], 0, 'g', 20);
  s << str3.arg(0.0, 0, 'g', 20)
            .arg(0.0, 0, 'g', 20)
            .arg(1/_slow[_layerCount-1], 0, 'g', 20) << flush;
}

void QGpCoreWave::RefractionDippingModel::toStream ( QDataStream & s ) const

QString QGpCoreWave::RefractionDippingModel::toString ( ) const [virtual]

Returns model as a string

Implements QGpCoreWave::GeophysicalModel.

References str, and TRACE.

{
  TRACE;
  QString str;
  str=QString("%1 layers, width %2 m\n").arg(_layerCount).arg(fabs(_xR-_xL));
  str+="  Left depth(m) Right depth(m)  Velocity(m/s)\n";
  for(int i=0;i<_layerCount-1;i++) {
    str+=QString::number(_depthL[i],'f',1).leftJustified(15,' ',true);
    str+=QString::number(_depthR[i],'f',1).leftJustified(15,' ',true);
    str+=QString::number(1.0/_slow[i],'f',1).leftJustified(15,' ',true);
    str+="\n";
  }
  str+=QString("-").leftJustified(15,' ',true);
  str+=QString("-").leftJustified(15,' ',true);
  str+=QString::number(1/_slow[_layerCount-1],'f',1).leftJustified(14,' ',true);
  str+="\n";
  return str;
}

double QGpCoreWave::RefractionDippingModel::travelTime	(	const TiltNode &	src,
		const TiltNode &	rec,
		int &	minTimeLayer
	)

Return travel time between two nodes src and rec. minTimeLayer is set to the deepest layer.

References TRACE.

Referenced by QGpCompatibility::CompatHodoChroneData::calculate(), QGpCoreWave::RefractionFactory::calculate(), RefraReader::parse(), and ray().

{
  TRACE;
  const TiltPath * srcPaths, * recPaths;
  minTimeLayer=0;
  if(!node2paths(src, rec, srcPaths, recPaths)) return 0.0;
  return travelTime(srcPaths, recPaths, minTimeLayer);
}

double QGpCoreWave::RefractionDippingModel::xLeft ( ) const [inline]

Referenced by QGpCoreWave::TiltNode::init().

{return _xL;}

double QGpCoreWave::RefractionDippingModel::xRight ( ) const [inline]

{return _xR;}

The documentation for this class was generated from the following files:

QGpCoreWave/RefractionDippingModel.h
QGpCoreWave/RefractionDippingModel.cpp

Classes

Public Types

Public Member Functions

Static Public Member Functions

Detailed Description

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation