Calculate Rayleigh ellipticity curves. More...
#include <Ellipticity.h>
Public Member Functions | |
| Ellipticity (int nModes, const QVector< double > *x) | |
| double | peakMisfit (int modeIndex, const RealStatisticalValue &val, Dispersion &disp, Rayleigh *waveModel) |
| QList< double > | peaks (int modeIndex, Dispersion &disp, Rayleigh *waveModel) |
Calculate Rayleigh ellipticity curves.
The main function used to calculate the ellipticity curve is calculate().
Experience has shown that a 1e-10 precision on dc is necessary to refine the ellipticity down to 1e-3 The precision used to compute the original curve must be the same as for the refines.
| QGpCoreWave::Ellipticity::Ellipticity | ( | int | nModes, |
| const QVector< double > * | x | ||
| ) | [inline] |
: ModalStorage(nModes, x) {}
| double QGpCoreWave::Ellipticity::peakMisfit | ( | int | modeIndex, |
| const RealStatisticalValue & | val, | ||
| Dispersion & | disp, | ||
| Rayleigh * | waveModel | ||
| ) |
Find the frequency of the closest peak to val and return the misfit
References QGpCoreTools::endl(), QGpCoreTools::Value< numberType >::isValid(), QGpCoreTools::StatisticalValue< numberType >::mean(), misfit(), QGpCoreWave::ModalStorage::mode(), QGpCoreWave::RayleighTemplate< RealType >::model(), QGpCoreWave::ModalStorage::refineSort(), QGpCoreWave::Seismic1DModel::roughFundamentalFrequency(), QGpCoreTools::StatisticalValue< numberType >::stddev(), QGpCoreTools::tr(), TRACE, QGpCoreWave::ModalStorage::x(), and QGpCoreWave::ModalStorage::xCount().
Referenced by QGpCoreWave::EllipticityFactory::peakMisfit().
{
TRACE;
double f0, devf0, misfit=1e99;
f0=val.mean()*2.0*M_PI;
if(val.stddev()>0) devf0=1.0/ (val.stddev()*2.0*M_PI); else devf0=1.0/f0;
RootSolver<Rayleigh> solver(waveModel);
RealValue * point=mode(modeIndex);
int n=xCount();
int i,j;
for(j=1; j<n && !point[j].isValid(); j++) {}
if(j==n) return 0.0;
for(i=j+1; i<n && !point[i].isValid(); i++ ) {}
if(i==n) return 0.0;
double lastSlope=fabs(point[i].value())-fabs(point[j].value());
j=i;
for(i++; i<n; i++ ) {
if(point[i].isValid()) {
double slope=fabs(point[i].value())-fabs(point[j].value());
j=i;
if(slope<=0 && lastSlope>0) {
double diff;
// Found a maximum in the sampled curve that is at distance > 0.01 the ellipticity of the half space
if(((diff=(f0-x(i))*devf0)>=0 && diff<misfit) ||
((diff=(x(i-2)-f0)*devf0)>=0 && diff<misfit) ||
(x(i-2)<f0 && f0<x(i))) {
// There's a chance to get a better misfit
bool badSampling;
double xr=refineMax(modeIndex, i, disp, solver, badSampling);
//App::stream() << tr("Refined peak frequency=%1").arg(0.5*xr/M_PI) << endl;
// The peak is sufficiently refined, adjust the misfit if necessary
diff=fabs(xr-f0)*devf0;
if(diff<misfit) misfit=diff;
if(badSampling)
App::stream() << tr("** Warning ** : probably a missing peak, bad initial sampling") << endl;
}
}
lastSlope=slope;
}
}
if(misfit==1e99) {
// monotonous curve, evaluate fundamental frequency by V0/4H
double roughFreq=2.0 * M_PI * waveModel->model()->roughFundamentalFrequency();
//App::stream() << tr("Rough fundamental frequency=%1").arg(0.5 * roughFreq/M_PI) << endl;
if(roughFreq>=x(0) && roughFreq<=x(n-1)) {
App::stream() << tr("** Warning ** : estimated rough fundamental frequency inside range") << endl;
double midRange=0.5*(x(0)+x(n-1));
if(roughFreq>midRange)
misfit=fabs(x(n-1)-f0)*devf0;
else
misfit=fabs(x(0)-f0)*devf0;
} else
misfit=fabs(roughFreq-f0)*devf0;
}
refineSort();
disp.refineSort();
return floor(misfit*1000+0.5)*0.001; // frequency is estimated down to 1e-3, avoid classification within precision limits
}
| QList< double > QGpCoreWave::Ellipticity::peaks | ( | int | modeIndex, |
| Dispersion & | disp, | ||
| Rayleigh * | waveModel | ||
| ) |
Return a list of the frequency of all peaks within the sampled range
References QGpCoreTools::endl(), QGpCoreTools::Value< numberType >::isValid(), QGpCoreWave::ModalStorage::mode(), QGpCoreTools::tr(), QGpCoreTools::Value< numberType >::value(), and QGpCoreWave::ModalStorage::xCount().
Referenced by EllipticityReader::parse().
{
QList<double> peaks;
double fac=0.5/M_PI;
RootSolver<Rayleigh> solver(waveModel);
RealValue * point=mode(modeIndex);
int n=xCount();
int i,j;
for(j=1; j<n && !point[j].isValid(); j++) {}
if(j==n) return peaks;
for(i=j+1; i<n && !point[i].isValid(); i++ ) {}
if(i==n) return peaks;
double lastSlope=point[i].value()-point[j].value();
j=i;
for(i++; i<n; i++ ) {
if(point[i].isValid()) {
double slope=fabs(point[i].value())-fabs(point[j].value());
j=i;
if(slope<=0 && lastSlope>0 && lastSlope-slope>1e-10) {
// Found a maximum in the sampled curve that is at distance > 0.01 the ellipticity of the half space
bool badSampling;
double xr=refineMax(modeIndex, i, disp, solver, badSampling);
peaks << fac*xr;
if(badSampling)
App::stream() << tr("** Warning ** : probably a missing peak, bad initial sampling") << endl;
}
lastSlope=slope;
}
}
return peaks;
}