All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines
Calculate Love and Rayleigh dispersion curves. More...
#include <Dispersion.h>
Public Member Functions | |
| bool | calculate (Rayleigh *waveModel, Ellipticity *ell=0, QAtomicInt *terminated=0) |
| bool | calculate (Love *waveModel, QAtomicInt *terminated=0) |
| QGPCOREWAVE_EXPORT ModalStorage * | deltaK () const |
| QGPCOREWAVE_EXPORT | Dispersion (int nModes, const QVector< double > *x) |
| QGPCOREWAVE_EXPORT bool | refine (int modeIndex, int lastIndex, double omega, RootSolver< Rayleigh > &modelS, Ellipticity &ell) |
| QGPCOREWAVE_EXPORT void | setGroupSlowness () |
| void | setPrecision (double p) |
Calculate Love and Rayleigh dispersion curves.
The main function used to calculate the dispersion curve is calculate().
| QGpCoreWave::Dispersion::Dispersion | ( | int | nModes, |
| const QVector< double > * | x | ||
| ) |
References TRACE.
: ModalStorage(nModes, x) { TRACE; _precision=1e-7; // Default precision without ellipticity }
| bool QGpCoreWave::Dispersion::calculate | ( | Rayleigh * | waveModel, |
| Ellipticity * | ell = 0, |
||
| QAtomicInt * | terminated = 0 |
||
| ) | [inline] |
Referenced by QGpCoreWave::DispersionFactory::calculate(), dispersion_curve_love_(), dispersion_curve_rayleigh_(), EllipticityReader::parse(), and DispersionReader::parse().
{
return calculate<double, Rayleigh>(waveModel, ell, terminated);
}
| bool QGpCoreWave::Dispersion::calculate | ( | Love * | waveModel, |
| QAtomicInt * | terminated = 0 |
||
| ) | [inline] |
{
return calculate<double, Love>(waveModel, 0, terminated);
}
| ModalStorage * QGpCoreWave::Dispersion::deltaK | ( | ) | const |
Returns a new ModalStorage with wavenumber differences between modes
References QGpCoreWave::ModalStorage::ModalStorage(), QGpCoreWave::ModalStorage::mode(), QGpCoreWave::ModalStorage::modeCount(), QGpCoreTools::Value< numberType >::setValid(), QGpCoreTools::Value< numberType >::setValue(), TRACE, QGpCoreWave::ModalStorage::x(), and QGpCoreWave::ModalStorage::xCount().
Referenced by DispersionReader::parse().
{
TRACE;
ModalStorage * s=new ModalStorage(modeCount()-1, x());
for(int im=modeCount()-1; im>0; im--) {
const RealValue * m1=mode(im);
const RealValue * m2=mode(im-1);
RealValue * dk=s->mode(im-1);
for(int is=xCount()-1; is>=0; is--) {
if(m1[is].isValid() && m2[is].isValid()) {
dk[is].setValue(x(is)*(m2[is].value()-m1[is].value()));
dk[is].setValid(true);
} else {
dk[is].setValid(false);
}
}
}
return s;
}
| bool QGpCoreWave::Dispersion::refine | ( | int | modeIndex, |
| int | lastIndex, | ||
| double | omega, | ||
| RootSolver< Rayleigh > & | solver, | ||
| Ellipticity & | ell | ||
| ) |
Used by the ellipticity computation, it allows the refinement of the ellipticity peak. The method for dispersion computation is the same as calculate(). It is just a resume to a particular frequency sample. lastIndex is the x index of the fixed and already computed curve. It avoids searching in the x list every time this function is called. lastIndex must be less than xCount()-1.
References QGpCoreWave::Seismic1DModel::checkVelocityInversion(), QGpCoreWave::RayleighTemplate< RealType >::ellipticity(), QGpCoreTools::endl(), QGpCoreTools::RootSolverTemplate< double, FunctionClass >::inversePolarity(), QGpCoreTools::RootSolverTemplate< double, FunctionClass >::lower(), QGpCoreWave::Seismic1DModel::maxSlowR(), QGpCoreWave::Seismic1DModel::minSlowS(), QGpCoreWave::ModalStorage::mode(), modeIndex, QGpCoreWave::RayleighTemplate< RealType >::model(), QGpCoreTools::RootSolverTemplate< double, FunctionClass >::neville(), QGpCoreTools::RootSolverTemplate< double, FunctionClass >::precision(), QGpCoreWave::ModalStorage::refineAdd(), SAFE_UNINITIALIZED, QGpCoreTools::RootSolverTemplate< double, FunctionClass >::searchDown(), QGpCoreTools::RootSolverTemplate< double, FunctionClass >::searchStep(), QGpCoreTools::RootSolverTemplate< double, FunctionClass >::searchUp(), QGpCoreWave::RayleighTemplate< RealType >::setOmega(), QGpCoreTools::RootSolverTemplate< double, FunctionClass >::setPolarity(), QGpCoreTools::RootSolverTemplate< double, FunctionClass >::setPrecision(), QGpCoreTools::RootSolverTemplate< double, FunctionClass >::setRelativeStep(), QGpCoreWave::ModalRefine::setValid(), QGpCoreWave::ModalRefine::setValue(), QGpCoreTools::tr(), TRACE, QGpCoreTools::RootSolverTemplate< double, FunctionClass >::upper(), QGpCoreTools::Value< numberType >::value(), QGpCoreTools::RootSolverTemplate< double, FunctionClass >::values(), QGpCoreWave::ModalStorage::x(), QGpCoreWave::ModalStorage::xCount(), and QGpCoreWave::RayleighTemplate< RealType >::y().
{
TRACE;
ASSERT(lastIndex<xCount()-1);
ModalRefine& dispRef=refineAdd(omega);
ModalRefine& ellRef=ell.refineAdd(omega);
RealValue * values=mode(0)+lastIndex;
Rayleigh * waveModel=solver.values();
const Seismic1DModel * model=waveModel->model();
double minSlow=model->minSlowS();
double maxSlow, curSlow, lastSupBound=model->maxSlowR(), newCurSlow;
SAFE_UNINITIALIZED(curSlow,0);
bool forbitVelocityInversion=model->checkVelocityInversion()==-1;
// Initialize equation's polarity
waveModel->setOmega(0.05); // the sign of polarity is the same at all frequencies
solver.setPolarity(model->maxSlowR());
for(int iMode=0; iMode <= modeIndex; iMode++ ) {
solver.setPrecision(_precision);
if(iMode==0) {
maxSlow=model->maxSlowR();
} else {
maxSlow=curSlow;
solver.inversePolarity();
}
while(true) {
waveModel->setOmega(omega);
//printf("\nFrequency=%.15lg Hz\n\n",0.5 * omega/M_PI);
bool errorDetected=false;
curSlow=values[iMode*xCount()+1].value();
if( ! solver.searchDown(curSlow, minSlow, maxSlow) ) {
if(iMode==0) {
// for fundamental mode this condition can never occur unless a mode jumping occured
App::stream() << tr( "** Warning ** : mode jumping for mode %1 (end), reducing "
"step ratio to %2" ).arg(iMode).arg(Number::toDouble(solver.searchStep() * 0.1) ) << endl;
errorDetected=true;
} else {
dispRef.setValid(iMode, false);
ellRef.setValid(iMode, false);
curSlow=minSlow;
}
} else {
solver.neville();
newCurSlow=solver.lower();
//printf("%lg %lg-->%lg %lg\n",minSlow, curSlow, newCurSlow, maxSlow);
waveModel->y(0.5 * (solver.upper() + newCurSlow) );
ellRef.setValue(iMode, waveModel->ellipticity());
if(newCurSlow - curSlow > _precision) {
if(forbitVelocityInversion) {
errorDetected=true;
App::stream() << tr( "** Warning ** : retrograde dispersion not allowed (mode %1), "
"reducing step ratio to %2" ).arg(iMode).arg(Number::toDouble(solver.searchStep() * 0.1) ) << endl;
} else {
// Inversions are allowed, check at last point that no mode exist at a higher slowness
waveModel->setOmega(x( lastIndex) );
solver.setRelativeStep(solver.searchStep() * 0.1);
if(solver.searchUp(lastSupBound, maxSlow) ) {
// recalculating curve because of mode jumping
App::stream() << tr( "** Warning ** : mode jumping for mode %1 (middle), "
"reducing step ratio to %2" ).arg(iMode).arg(Number::toDouble(solver.searchStep() * 0.1) ) << endl;
errorDetected=true;
}
solver.setRelativeStep(solver.searchStep() * 10.0);
}
}
if(!errorDetected) {
lastSupBound=solver.upper();
curSlow=newCurSlow;
dispRef.setValue(iMode, curSlow);
// When calculating the dispersion curve, stepRatio is positive to search downwards
// (From maxSlow towards minSlow)
// Here we want to search carefully upwards: from curSlow to maxSlow
// If at least one root is found, this proves that mode jumping occured
if(curSlow > minSlow) {
// The curSlow was calculated using _solve and _x2 is the value returned by
// Solve (slightly greater than the true root).
// Here we need the _x1 value, e.g. slightly less than the true root, in
// order to be sure that the root eventually found will not be the same
// as curSlow
curSlow=solver.upper();
}
// Contrary to normal dispersion curve computation, here we know the next point on the curve
// Hence, we simulate its computation and we must find the same number, if not mode jumping occured.
waveModel->setOmega(x( lastIndex) );
//printf("Check refine jumpig: frequency=%.15lg Hz\n",0.5 * x(lastIndex)/M_PI);
double correctSlow=values[iMode*xCount()].value();
bool ok;
if(iMode==0) {
ok=solver.searchDown(curSlow, minSlow, maxSlow);
} else {
double lastMaxSlow=values[(iMode-1)*xCount()].value();
if(curSlow>lastMaxSlow) curSlow=lastMaxSlow;
ok=solver.searchDown(curSlow, minSlow, lastMaxSlow);
}
if(ok) {
solver.neville();
newCurSlow=solver.lower();
}
//printf("%lg %lg-->%lg %lg=%lg %lg %lg\n",minSlow, curSlow, newCurSlow, lastMaxSlow, correctSlow, newCurSlow-correctSlow, _precision);
if( !ok || fabs(newCurSlow-correctSlow)>_precision) {
ok=true;
App::stream() << tr( "** Warning ** : mode jumping for mode %1 (refine), reducing "
"step ratio to %2" ).arg(iMode).arg(Number::toDouble(solver.searchStep() * 0.1) ) << endl;
errorDetected=true;
}
}
}
if(errorDetected) {
// recalculating curve because of mode jumping
if(solver.searchStep() < 1e-4) {
App::stream() << tr( "** Warning ** : dispersion curve impossible to obtain "
"for mode %1" ).arg(iMode) << endl;
return false;
}
solver.setPrecision(solver.precision() * 0.1);
solver.setRelativeStep(solver.searchStep() * 0.1);
} else break;
}
}
return true;
}
Convert phase slownesses into group slownesses
References QGpCoreWave::ModalStorage::mode(), QGpCoreWave::ModalStorage::modeCount(), QGpCoreTools::Value< numberType >::setValid(), TRACE, QGpCoreTools::Value< numberType >::value(), QGpCoreWave::ModalStorage::x(), and QGpCoreWave::ModalStorage::xCount().
Referenced by QGpCoreWave::DispersionFactory::calculate(), dispersion_curve_love_(), dispersion_curve_rayleigh_(), and DispersionReader::parse().
{
TRACE;
int nx=xCount()-1;
RealValue tmp[ xCount() ];
for(int im=0;im<modeCount();im++) {
RealValue * point=mode(im);
for(int i=1;i<nx;i++) {
if(point[i-1].isValid() && point[i].isValid() && point[i+1].isValid()) {
double val=point[i].value();
double derslow=0.5*((val-point[i-1].value())/(x(i)-x(i-1))+
(point[i+1].value()-val)/(x(i+1)-x(i)));
derslow*=x(i);
tmp[ i ]=val+derslow;
} else {
tmp[ i ].setValid(false);
}
}
point[0].setValid(false);
for(int i=1;i<nx;i++) {
point[i]=tmp[i];
}
point[nx].setValid(false);
}
}
| void QGpCoreWave::Dispersion::setPrecision | ( | double | p | ) | [inline] |
Set precision in slowness for dispersion computation. If you need more than 1e-16, you must use the big number version of calculate().
Referenced by EllipticityReader::parse().
{_precision=p;}
