All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines
Implementation of Rayleigh surface waves. More...
#include <Rayleigh.h>
Public Member Functions | |
| double | ellipticity () |
| const Seismic1DModel * | model () const |
| RayleighTemplate (const Seismic1DModel *model) | |
| void | setOmega (double o) |
| RealType | y (RealType slowness) |
Implementation of Rayleigh surface waves.
It is a layered model able to calculate the roots corresponding to the theoretical Rayleigh dispersion curve.
| QGpCoreWave::RayleighTemplate< RealType >::RayleighTemplate | ( | const Seismic1DModel * | model | ) | [inline] |
References QGpCoreWave::RayleighTemplate< RealType >::model().
{_model=model;}
| double QGpCoreWave::RayleighTemplate< RealType >::ellipticity | ( | ) | [inline] |
Return the ellipticity for current frequency (omega). Once y(slowness) has been run through a solve() call, you can call this function to get the ellipticity, ratio of eigenfunctions for vertical and horizontal components.
Referenced by QGpCoreWave::Dispersion::refine().
{return Number::toDouble(iF1214/F1213);}
| const Seismic1DModel* QGpCoreWave::RayleighTemplate< RealType >::model | ( | ) | const [inline] |
Referenced by QGpCoreWave::Ellipticity::peakMisfit(), QGpCoreWave::RayleighTemplate< RealType >::RayleighTemplate(), and QGpCoreWave::Dispersion::refine().
{return _model;}
| void QGpCoreWave::RayleighTemplate< RealType >::setOmega | ( | double | o | ) | [inline] |
Referenced by QGpCoreWave::ModalCurve::determinantMisfit(), DispersionReader::parse(), QGpCoreWave::Dispersion::refine(), and QGpCoreWave::RayleighFunction::value().
{_omega=o;}
| RealType QGpCoreWave::RayleighTemplate< RealType >::y | ( | RealType | slowness | ) |
Calculates the Rayleigh surface wave equation. This function is automatically called by the Neville Solver. Implementation of Dunkin's method (BSSA 1965), almost the same as Herrmann. Modification of the formulation by the author of this C++ class to optimize the CPU time needed.
Normally called only by the root solver.
References QGpCoreTools::cos(), QGpCoreTools::exp(), G1313, G1334, G1414, G2314, G2424, G2434, G3413, G3424, G3434, iG1424, iG1434, iG2414, iG3414, QGpCoreTools::sin(), and QGpCoreTools::sqrt().
Referenced by QGpCoreWave::ModalCurve::determinantMisfit(), DispersionReader::parse(), QGpCoreWave::Dispersion::refine(), and QGpCoreWave::RayleighFunction::value().
{
RealType tmpAbs;
//printf("f=%.20lg s=%.20lg v=%.20lg\n",_omega/(2*M_PI),slowness,1.0/slowness);
RealType w2=(RealType)_omega*(RealType)_omega;
RealType k=(RealType)_omega*slowness;
RealType k2=k*k;
RealType inv_k2=(RealType)1.0/k2;
// Half space
register int i=_model->layerCount()-1;
RealType ka=(RealType)_omega*(RealType)_model->slowP(i);
RealType hn2=k2-ka*ka;
tmpAbs=hn2;
if(tmpAbs<0) tmpAbs=-tmpAbs;
RealType hn=sqrt(tmpAbs);
RealType kb=(RealType)_omega*(RealType)_model->slowS(i);
RealType kn2=k2-kb*kb;
tmpAbs=kn2;
if(tmpAbs<0) tmpAbs=-tmpAbs;
RealType kn=sqrt(tmpAbs);
RealType ln=k2+kn2;
RealType mu=(RealType)_model->mu(i);
RealType c1=hn*kn;
// Factor used to balance the values in the plane (k,s)
//RealType invFac=w2*k2*slowness;
RealType invFac=w2;
RealType fac=1.0/invFac;
RealType F1212=mu*fac*(ln*ln-4.0*k2*c1);
F1213=hn*(ln-2.0*k2);
iF1214=k*(ln-2.0*c1);
// iF1223=iF1214
RealType F1224=kn*(2.0*k2-ln);
RealType F1234=(k2-c1)/mu;
for(i--;i>=0;i--) {
ka=(RealType)_omega*(RealType)_model->slowP(i);
hn2=k2-ka*ka;
tmpAbs=hn2;
if(tmpAbs<0) tmpAbs=-tmpAbs;
hn=sqrt(tmpAbs);
hn2*=inv_k2;
kb=(RealType)_omega*(RealType)_model->slowS(i);
kn2=k2-kb*kb;
tmpAbs=kn2;
if(tmpAbs<0) tmpAbs=-tmpAbs;
kn=sqrt(tmpAbs);
kn2*=inv_k2;
// Compute the hyperbolic sin and cos
RealType SH, CH, exphn;
//RealType he;
RealType dn=(RealType)_model->h(i);
exphn=hn*dn;
if(hn2>0.0) {
//he=exphn;
RealType exp2;
//if(exphn<21.2) { // See Love for explanations
if(exphn<115) {
exphn=exp(-exphn);
exp2=exphn*exphn;
} else {
exphn=0.0;
exp2=0.0;
}
SH=0.5*k/hn*(1.0-exp2);
CH=0.5*(1.0+exp2);
} else if(hn2<0.0) {
//he=0.0;
// sinh is imaginary but hn as well, thus SH is real
SH=k/hn*sin(exphn);
CH=cos(exphn);
exphn=1.0;
} else {
//he=0.0;
SH=0.0;
CH=1.0;
exphn=1.0;
}
RealType SK, CK, expkn;
//RealType ke;
expkn=kn*dn;
if(kn2>0.0) {
//ke=expkn;
RealType exp2;
if(expkn<21.2) { // See Love for explanations (double precision)
//if(expkn<115) { // See Love for explanations (big number precise until 1e-50)
expkn=exp(-expkn);
exp2=expkn*expkn;
} else {
expkn=0.0;
exp2=0.0;
}
SK=0.5*k/kn*(1.0-exp2);
CK=0.5*(1.0+exp2);
} else if(kn2<0.0) {
//ke=0.0;
// sinh is imaginary but kn as well, thus SK is real
SK=k/kn*sin(expkn);
CK=cos(expkn);
expkn=1;
} else {
//ke=0.0;
SK=0.0;
CK=1.0;
expkn=1.0;
}
// In G, all terms that does not contains either CHCK,... we must multiply by expCorr
RealType expCorr=exphn*expkn;
// Each cosh or sinh has to be multiplied by either exp(he) or exp(ke)
// he is always > ke, thus exp(he)>exp(ke).
// In each term, CHCK, CHSK, ... there is a factor exp(he+ke)
RealType CHCK=CH*CK;
RealType SHSK=SH*SK;
RealType CHSK=CH*SK;
RealType SHCK=SH*CK;
// CHCK is always real
// Other may be imaginary or complex
// Intermediate variable, avoid dupplicate operations
RealType gn=2.0*k2/(kb*kb);
RealType gn2=gn*gn;
RealType adim1=gn2-2.0*gn+1.0;
RealType adim2=hn2*kn2;
RealType adim3=gn2+adim1;
RealType adim4=1.0-gn;
RealType adim5=gn2*adim2;
RealType CHCK1=expCorr-CHCK;
RealType CHCK2=2.0*CHCK1;
RealType SHCK1=gn*hn2*SHCK;
c1=(RealType)_model->rho(i)*w2/k;
RealType c2=1.0/c1;
// Second order subdeterminants of propagator matrix
// i means that the subdeterminant is imaginary
RealType G1212=adim3*CHCK-(adim1+adim5)*SHSK-(adim3-1)*expCorr;
RealType G1213=c2*(CHSK-hn2*SHCK);
RealType iG1214=c2*((adim1-gn2)*CHCK1+(adim4-gn*adim2)*SHSK);
//#define iG1223 iG1214 --> R1223
RealType G1224=c2*(kn2*CHSK-SHCK);
RealType G1234=c2*c2*(CHCK2+(1+adim2)*SHSK);
RealType G1312=c1*(gn2*kn2*CHSK-adim1*SHCK);
#define G1313 CHCK
RealType iG1314=adim4*SHCK+gn*kn2*CHSK;
//#define iG1323 iG1314 --> R1223
RealType G1324=kn2*SHSK;
#define G1334 G1224
RealType iG1412=c1*((adim1-adim4)*(adim4-gn)*CHCK1+(adim4*adim1-gn*adim5)*SHSK);
RealType iG1413=SHCK1+adim4*CHSK;
RealType G1423=CHCK-G1212;
#define G1414 expCorr+G1423
#define iG1424 iG1314
#define iG1434 iG1214
//#define iG2312 iG1412
//#define iG2313 iG1413
#define G2314 G1423
//#define G2323 G1414 --> R1223
//#define iG2324 iG1314
//#define iG2334 iG1214
RealType G2412= c1*(adim1*CHSK-gn*SHCK1);
RealType G2413=hn2*SHSK;
#define iG2414 iG1413
//#define iG2423 iG1413 --> R1223
#define G2424 G1313
#define G2434 G1213
RealType G3412=c1*c1*(gn2*adim1*CHCK2+(adim1*adim1+gn2*adim5)*SHSK);
#define G3413 G2412
#define iG3414 iG1412
//#define iG3423 iG1412 --> R1223
#define G3424 G1312
#define G3434 G1212
/*
Generic formula:
R12cd=F1212*G12cd+F1213*G13cd+iF1214*iG14cd+iF1223*iG23cd+F1224*G24cd+F1234*G34cd
*/
RealType R1212=F1212*G1212+fac*(F1213*G1312-2.0*iF1214*iG1412+F1224*G2412-F1234*G3412);
RealType R1213=invFac*F1212*G1213+F1213*G1313+2.0*iF1214*iG1413-F1224*G2413+F1234*G3413;
RealType iR1214= invFac*F1212*iG1214+F1213*iG1314+iF1214*(G1414+G2314)-F1224*iG2414+F1234*iG3414;
// R1223=invFac*F1212*iG1223+F1213*iG1323+iF1214*(G1423+G2323)-F1224*iG2423+F1234*iG3423;
RealType R1224=invFac*F1212*G1224-F1213*G1324-2.0*iF1214*iG1424+F1224*G2424+F1234*G3424;
RealType R1234=F1213*G1334-invFac*F1212*G1234-2.0*iF1214*iG1434+F1224*G2434+F1234*G3434;
// Normalize R before copying to F
RealType maxR=R1212;
if(maxR<0) maxR=-maxR;
tmpAbs=R1213;
if(tmpAbs<0) tmpAbs=-tmpAbs;
if(tmpAbs>maxR) maxR=tmpAbs;
tmpAbs=iR1214;
if(tmpAbs<0) tmpAbs=-tmpAbs;
if(tmpAbs>maxR) maxR=tmpAbs;
tmpAbs=R1224;
if(tmpAbs<0) tmpAbs=-tmpAbs;
if(tmpAbs>maxR) maxR=tmpAbs;
tmpAbs=R1234;
if(tmpAbs<0) tmpAbs=-tmpAbs;
if(tmpAbs>maxR) maxR=tmpAbs;
if(maxR>1.0e5) {
maxR=1.0e5/maxR;
F1212=R1212*maxR;
F1213=R1213*maxR;
iF1214=iR1214*maxR;
F1224=R1224*maxR;
F1234=R1234*maxR;
}
else {
F1212=R1212;
F1213=R1213;
iF1214=iR1214;
F1224=R1224;
F1234=R1234;
}
}
return F1212;
}
