discur3frack.h

弹性波中瑞雷波的正演计算。求解多层瑞雷波频散曲线。计算速度快。（专业小软件）

//////////////////////////fadhrtf()三层中间裂缝夹层计算///////////////////////////////////
#include"math.h" 
#include"jzz.h"
//本模型取第一层和第三层介质弹性参数一致
//vr瑞雷波速度
//vs第一层和第三层横波速度
//vp第一层和第三层纵波速度
//vf第二层介质波速
//den1为一、三层的介质密度，den2为中间层的介质密度
//h1,h2分别为第一层和第二层的介质厚度
//f源频率
double fadhrtf(double vr,double vs,double vp,double den1,double den2,double vf,double h1,double h2,double f)
{
	double hls[25],det,B,r,rp,rs,rf,Q,ep1,ep2,es1,es2,ef,A;	
	r=2.0*vs*vs/(vr*vr);
	rp=sqrt(1.0-vr*vr/(vp*vp));
	rs=sqrt(1.0-vr*vr/(vs*vs));
	rf=sqrt(vr*vr/(vf*vf)-1.0);
	Q=den2/den1;
	ep1=exp(2.0*3.1415926*f*rp*h1/vr);
	ep2=1.0/ep1;
	es1=exp(2.0*3.1415926*f*rs*h1/vr);
	es2=1.0/es1;
	ef=2.0*3.1415926*f*rf*h2/vr;
	A=rp*Q*cos(ef)-(pow(r-1,2)-r*r*rp*rs)*rf*sin(ef);
	B=rp*Q*sin(ef)+(pow(r-1,2)-r*r*rp*rs)*rf*cos(ef);
    hls[0]=r-1; hls[1]=r-1; hls[2]=-1.0*r*rs; hls[3]=r*rs; hls[4]=0;
	hls[5]=-1.0*r*rp; hls[6]=r*rp; hls[7]=r-1; hls[8]=r-1; hls[9]=0;
    hls[10]=rp*ep1;hls[11]=-1.0*rp*ep2;hls[12]=-1.0*es1;hls[13]=-1.0*es2;hls[14]=rf*A;
	hls[15]=(r-1)*ep1; hls[16]=(r-1)*ep2; hls[17]=-1.0*r*rs*es1; hls[18]=r*rs*es2; hls[19]=Q*B;
	hls[20]=-1.0*r*rp*ep1; hls[21]=r*rp*ep2; hls[22]=(r-1)*es1; hls[23]=(r-1)*es2; hls[24]=0.0;
    det=bsdet(hls,5);
	return(det);
   
}
		
//////////三层中间裂缝介质二分法计算/////////////
//fadhrt()三层固体介质瑞雷波频散方程二分法求解算法,求解需用到函数fadhrtf()
//a,b为根的求解区间，a<b，一般a大于第二层介质的横波速度，b小于一、三横波速度;
//h为二分法求解步长;
//eps控制精度;
//x一维数组，长度为m[0]，返回在区间[a,b]中搜度到的瑞雷波速度的解;
//m一维数组，长度为2，一般取m[0]=10,取m[1]=0;
//vs第一层和第三层横波速度
//vp第一层和第三层纵波速度
//vf第二层介质波速
//den1为一、三层的介质密度，den2为中间层的介质密度
//h1,h2分别为第一层和第二层的介质厚度
//f源频率		
		void fadhrt(double a1,double b,double h,double eps,double x[],int m[],
		double vs,double vp,double den1,double den2,double vf,double h1,double h2,double f)
		{ 
			int n,js;
			double z,y,z1,y1,z0,y0;
			n=0;
			z=a1;
			y=fadhrtf(z,vs,vp,den1,den2,vf,h1,h2,f);
			while((z<=b+h/2.0)&&(n!=m[0]))
			{
				if(fabs(y)<eps)
				{  
					n=n+1;x[n-1]=z;
					z=z+h/2.0;y=fadhrtf(z,vs,vp,den1,den2,vf,h1,h2,f);
				}
				
				else
				{  
					z1=z+h;y1=fadhrtf(z1,vs,vp,den1,den2,vf,h1,h2,f);      
					if(fabs(y1)<eps)
					{n=n+1;x[n-1]=z1;
					z=z1+h/2.0;y=fadhrtf(z,vs,vp,den1,den2,vf,h1,h2,f);
					}
					else if(y*y1>0.0)
					{
						y=y1;z=z1;}
					else
					{
						js=0;
						while(js==0)
						{
							if(fabs(z1-z)<eps)
							{
								n=n+1;x[n-1]=(z1+z)/2.0;
								z=z1+h/2.0;y=fadhrtf(z,vs,vp,den1,den2,vf,h1,h2,f);
								js=1;
							}
							else
							{
								z0=(z1+z)/2.0;y0=fadhrtf(z0,vs,vp,den1,den2,vf,h1,h2,f);
								if(fabs(y0)<eps)
								{x[n]=z0;n=n+1;js=1;
								z=z0+h/2.0;
								y=fadhrtf(z,vs,vp,den1,den2,vf,h1,h2,f);
								}
								else if ((y*y0)<0.0)
								{z=z0;y=y0;} 
								else {z=z0;y=y0;}
							}
						}
					}
				}
			}
			m[1]=n;
			return;
		}