00001 /*************************************************************************** 00002 ** 00003 ** This file is part of QGpCoreWave. 00004 ** 00005 ** This library is free software; you can redistribute it and/or 00006 ** modify it under the terms of the GNU Lesser General Public 00007 ** License as published by the Free Software Foundation; either 00008 ** version 2.1 of the License, or (at your option) any later version. 00009 ** 00010 ** This file is distributed in the hope that it will be useful, but WITHOUT 00011 ** ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 ** FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 00013 ** License for more details. 00014 ** 00015 ** You should have received a copy of the GNU Lesser General Public 00016 ** License along with this library; if not, write to the Free Software 00017 ** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00018 ** 00019 ** See http://www.geopsy.org for more information. 00020 ** 00021 ** Created : 2004-04-28 00022 ** Authors : 00023 ** Marc Wathelet 00024 ** Marc Wathelet (ULg, Liège, Belgium) 00025 ** Marc Wathelet (LGIT, Grenoble, France) 00026 ** 00027 ***************************************************************************/ 00028 00029 #ifndef RAYLEIGH_H 00030 #define RAYLEIGH_H 00031 00032 #include <QGpCoreTools.h> 00033 00034 #include "QGpCoreWaveDLLExport.h" 00035 #include "Seismic1DModel.h" 00036 00037 namespace QGpCoreWave { 00038 00039 template<class RealType> 00040 class RayleighTemplate 00041 { 00042 public: 00043 RayleighTemplate(const Seismic1DModel * model) {_model=model;} 00044 00045 const Seismic1DModel * model() const {return _model;} 00046 00047 void setOmega(double o) {_omega=o;} 00048 double ellipticity() {return Number::toDouble(iF1214/F1213);} 00049 RealType y(RealType slowness); 00050 private: 00051 const Seismic1DModel * _model; 00052 /* Dunkin determinants for the surperficial layer (once y(double) has run, 00053 else undetermined) to compute ellipticity 00054 */ 00055 RealType F1213; 00056 RealType iF1214; 00057 double _omega; 00058 }; 00059 00060 // Usual precision for efficient computations 00061 typedef RayleighTemplate<double> Rayleigh; 00062 #ifdef BIGNUM 00063 // Arbitrary precision for special cases or experimental purposes 00064 typedef RayleighTemplate<mp_real> RayleighBigNum; 00065 #endif 00066 00067 template<class RealType> 00068 RealType RayleighTemplate<RealType>::y(RealType slowness) 00069 { 00070 RealType tmpAbs; 00071 //printf("f=%.20lg s=%.20lg v=%.20lg\n",_omega/(2*M_PI),slowness,1.0/slowness); 00072 RealType w2=(RealType)_omega*(RealType)_omega; 00073 RealType k=(RealType)_omega*slowness; 00074 RealType k2=k*k; 00075 RealType inv_k2=(RealType)1.0/k2; 00076 // Half space 00077 register int i=_model->layerCount()-1; 00078 RealType ka=(RealType)_omega*(RealType)_model->slowP(i); 00079 RealType hn2=k2-ka*ka; 00080 tmpAbs=hn2; 00081 if(tmpAbs<0) tmpAbs=-tmpAbs; 00082 RealType hn=sqrt(tmpAbs); 00083 RealType kb=(RealType)_omega*(RealType)_model->slowS(i); 00084 RealType kn2=k2-kb*kb; 00085 tmpAbs=kn2; 00086 if(tmpAbs<0) tmpAbs=-tmpAbs; 00087 RealType kn=sqrt(tmpAbs); 00088 RealType ln=k2+kn2; 00089 RealType mu=(RealType)_model->mu(i); 00090 RealType c1=hn*kn; 00091 // Factor used to balance the values in the plane (k,s) 00092 //RealType invFac=w2*k2*slowness; 00093 RealType invFac=w2; 00094 RealType fac=1.0/invFac; 00095 RealType F1212=mu*fac*(ln*ln-4.0*k2*c1); 00096 F1213=hn*(ln-2.0*k2); 00097 iF1214=k*(ln-2.0*c1); 00098 // iF1223=iF1214 00099 RealType F1224=kn*(2.0*k2-ln); 00100 RealType F1234=(k2-c1)/mu; 00101 for(i--;i>=0;i--) { 00102 ka=(RealType)_omega*(RealType)_model->slowP(i); 00103 hn2=k2-ka*ka; 00104 tmpAbs=hn2; 00105 if(tmpAbs<0) tmpAbs=-tmpAbs; 00106 hn=sqrt(tmpAbs); 00107 hn2*=inv_k2; 00108 kb=(RealType)_omega*(RealType)_model->slowS(i); 00109 kn2=k2-kb*kb; 00110 tmpAbs=kn2; 00111 if(tmpAbs<0) tmpAbs=-tmpAbs; 00112 kn=sqrt(tmpAbs); 00113 kn2*=inv_k2; 00114 // Compute the hyperbolic sin and cos 00115 RealType SH, CH, exphn; 00116 //RealType he; 00117 RealType dn=(RealType)_model->h(i); 00118 exphn=hn*dn; 00119 if(hn2>0.0) { 00120 //he=exphn; 00121 RealType exp2; 00122 //if(exphn<21.2) { // See Love for explanations 00123 if(exphn<115) { 00124 exphn=exp(-exphn); 00125 exp2=exphn*exphn; 00126 } else { 00127 exphn=0.0; 00128 exp2=0.0; 00129 } 00130 SH=0.5*k/hn*(1.0-exp2); 00131 CH=0.5*(1.0+exp2); 00132 } else if(hn2<0.0) { 00133 //he=0.0; 00134 // sinh is imaginary but hn as well, thus SH is real 00135 SH=k/hn*sin(exphn); 00136 CH=cos(exphn); 00137 exphn=1.0; 00138 } else { 00139 //he=0.0; 00140 SH=0.0; 00141 CH=1.0; 00142 exphn=1.0; 00143 } 00144 RealType SK, CK, expkn; 00145 //RealType ke; 00146 expkn=kn*dn; 00147 if(kn2>0.0) { 00148 //ke=expkn; 00149 RealType exp2; 00150 if(expkn<21.2) { // See Love for explanations (double precision) 00151 //if(expkn<115) { // See Love for explanations (big number precise until 1e-50) 00152 expkn=exp(-expkn); 00153 exp2=expkn*expkn; 00154 } else { 00155 expkn=0.0; 00156 exp2=0.0; 00157 } 00158 SK=0.5*k/kn*(1.0-exp2); 00159 CK=0.5*(1.0+exp2); 00160 } else if(kn2<0.0) { 00161 //ke=0.0; 00162 // sinh is imaginary but kn as well, thus SK is real 00163 SK=k/kn*sin(expkn); 00164 CK=cos(expkn); 00165 expkn=1; 00166 } else { 00167 //ke=0.0; 00168 SK=0.0; 00169 CK=1.0; 00170 expkn=1.0; 00171 } 00172 // In G, all terms that does not contains either CHCK,... we must multiply by expCorr 00173 RealType expCorr=exphn*expkn; 00174 // Each cosh or sinh has to be multiplied by either exp(he) or exp(ke) 00175 // he is always > ke, thus exp(he)>exp(ke). 00176 // In each term, CHCK, CHSK, ... there is a factor exp(he+ke) 00177 RealType CHCK=CH*CK; 00178 RealType SHSK=SH*SK; 00179 RealType CHSK=CH*SK; 00180 RealType SHCK=SH*CK; 00181 // CHCK is always real 00182 // Other may be imaginary or complex 00183 // Intermediate variable, avoid dupplicate operations 00184 RealType gn=2.0*k2/(kb*kb); 00185 RealType gn2=gn*gn; 00186 RealType adim1=gn2-2.0*gn+1.0; 00187 RealType adim2=hn2*kn2; 00188 RealType adim3=gn2+adim1; 00189 RealType adim4=1.0-gn; 00190 RealType adim5=gn2*adim2; 00191 RealType CHCK1=expCorr-CHCK; 00192 RealType CHCK2=2.0*CHCK1; 00193 RealType SHCK1=gn*hn2*SHCK; 00194 c1=(RealType)_model->rho(i)*w2/k; 00195 RealType c2=1.0/c1; 00196 00197 // Second order subdeterminants of propagator matrix 00198 // i means that the subdeterminant is imaginary 00199 RealType G1212=adim3*CHCK-(adim1+adim5)*SHSK-(adim3-1)*expCorr; 00200 RealType G1213=c2*(CHSK-hn2*SHCK); 00201 RealType iG1214=c2*((adim1-gn2)*CHCK1+(adim4-gn*adim2)*SHSK); 00202 //#define iG1223 iG1214 --> R1223 00203 RealType G1224=c2*(kn2*CHSK-SHCK); 00204 RealType G1234=c2*c2*(CHCK2+(1+adim2)*SHSK); 00205 RealType G1312=c1*(gn2*kn2*CHSK-adim1*SHCK); 00206 #define G1313 CHCK 00207 RealType iG1314=adim4*SHCK+gn*kn2*CHSK; 00208 //#define iG1323 iG1314 --> R1223 00209 RealType G1324=kn2*SHSK; 00210 #define G1334 G1224 00211 RealType iG1412=c1*((adim1-adim4)*(adim4-gn)*CHCK1+(adim4*adim1-gn*adim5)*SHSK); 00212 RealType iG1413=SHCK1+adim4*CHSK; 00213 RealType G1423=CHCK-G1212; 00214 #define G1414 expCorr+G1423 00215 #define iG1424 iG1314 00216 #define iG1434 iG1214 00217 //#define iG2312 iG1412 00218 //#define iG2313 iG1413 00219 #define G2314 G1423 00220 //#define G2323 G1414 --> R1223 00221 //#define iG2324 iG1314 00222 //#define iG2334 iG1214 00223 RealType G2412= c1*(adim1*CHSK-gn*SHCK1); 00224 RealType G2413=hn2*SHSK; 00225 #define iG2414 iG1413 00226 //#define iG2423 iG1413 --> R1223 00227 #define G2424 G1313 00228 #define G2434 G1213 00229 RealType G3412=c1*c1*(gn2*adim1*CHCK2+(adim1*adim1+gn2*adim5)*SHSK); 00230 #define G3413 G2412 00231 #define iG3414 iG1412 00232 //#define iG3423 iG1412 --> R1223 00233 #define G3424 G1312 00234 #define G3434 G1212 00235 00236 /* 00237 Generic formula: 00238 R12cd=F1212*G12cd+F1213*G13cd+iF1214*iG14cd+iF1223*iG23cd+F1224*G24cd+F1234*G34cd 00239 */ 00240 RealType R1212=F1212*G1212+fac*(F1213*G1312-2.0*iF1214*iG1412+F1224*G2412-F1234*G3412); 00241 RealType R1213=invFac*F1212*G1213+F1213*G1313+2.0*iF1214*iG1413-F1224*G2413+F1234*G3413; 00242 RealType iR1214= invFac*F1212*iG1214+F1213*iG1314+iF1214*(G1414+G2314)-F1224*iG2414+F1234*iG3414; 00243 // R1223=invFac*F1212*iG1223+F1213*iG1323+iF1214*(G1423+G2323)-F1224*iG2423+F1234*iG3423; 00244 RealType R1224=invFac*F1212*G1224-F1213*G1324-2.0*iF1214*iG1424+F1224*G2424+F1234*G3424; 00245 RealType R1234=F1213*G1334-invFac*F1212*G1234-2.0*iF1214*iG1434+F1224*G2434+F1234*G3434; 00246 00247 // Normalize R before copying to F 00248 RealType maxR=R1212; 00249 if(maxR<0) maxR=-maxR; 00250 tmpAbs=R1213; 00251 if(tmpAbs<0) tmpAbs=-tmpAbs; 00252 if(tmpAbs>maxR) maxR=tmpAbs; 00253 tmpAbs=iR1214; 00254 if(tmpAbs<0) tmpAbs=-tmpAbs; 00255 if(tmpAbs>maxR) maxR=tmpAbs; 00256 tmpAbs=R1224; 00257 if(tmpAbs<0) tmpAbs=-tmpAbs; 00258 if(tmpAbs>maxR) maxR=tmpAbs; 00259 tmpAbs=R1234; 00260 if(tmpAbs<0) tmpAbs=-tmpAbs; 00261 if(tmpAbs>maxR) maxR=tmpAbs; 00262 if(maxR>1.0e5) { 00263 maxR=1.0e5/maxR; 00264 F1212=R1212*maxR; 00265 F1213=R1213*maxR; 00266 iF1214=iR1214*maxR; 00267 F1224=R1224*maxR; 00268 F1234=R1234*maxR; 00269 } 00270 else { 00271 F1212=R1212; 00272 F1213=R1213; 00273 iF1214=iR1214; 00274 F1224=R1224; 00275 F1234=R1234; 00276 } 00277 } 00278 return F1212; 00279 } 00280 00281 /* 00282 kb, kn2, kn, 00283 From Herrmann (srfdis.f), for water layer on top 00284 if(llw.ne.1) then 00285 c----- 00286 c include water layer. 00287 c----- 00288 xka=omega/dble(a(1)) 00289 wvnop=wvno+xka 00290 wvnom=dabs(wvno-xka) 00291 ra=dsqrt(wvnop*wvnom) 00292 dpth=dble(d(1)) 00293 rho1=dble(rho(1)) 00294 p=ra*dpth 00295 beta=dble(b(1)) 00296 znul=1.0d-05 00297 call var(p,znul,ra,znul,wvno,xka,znul,dpth,w,cosp,exa) 00298 w0=-rho1*w 00299 dltar4=cosp*e(1) + w0*e(2) 00300 else 00301 dltar4=e(1) 00302 endif 00303 */ 00304 00305 } // namespace QGpCoreWave 00306 00307 #endif // RAYLEIGH_H
