discu2liqmode.h

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

#include"math.h"
// fadhrtf2liq()两层液体介质频散方程算法实现////////////////////////////
//vr瑞雷波速度，H第一层液体介质厚度，vp1第一层液体纵波速度，vp2第二层液体纵波速度
//den1第一层液体介质密度，den2第二层液体介质密度，f源频率
double fadhrtf2liq(double vr,double H,double vp1,double vp2,double den1,double den2,double f)
{
	///////////////可调参数/////////////////////////
	
	double out;	
	out=tan((2.0/vr)*3.1415926*f*H*sqrt(vr*vr/(vp1*vp1)-1.0))+(den2/den1)*sqrt(vr*vr/(vp1*vp1)-1.0)/sqrt(1.0-vr*vr/(vp2*vp2));			
 		
	return(out);
	
	
}
// fadhrt2liq()两层液体介质频散方程二分法求解算法,求解需用到函数fadhrtf2liq()
//a,b为根的求解区间，a<b，因为待求解的瑞雷波速度理论上介于第一层液体纵波速度和第二层液体纵波速度之间，所以一般取a>a1,b<a2;
//h为二分法求解步长;
//eps控制精度;
//x一维数组，长度为m[0]，返回在区间[a,b]中搜度到的瑞雷波速度的解;
//m一维数组，长度为2，一般取m[0]=10,取m[1]=0;
//H第一层液体介质厚度;
//a1第一层液体纵波速度，a2第二层液体纵波速度，一般取a1<a2;
//den1第一层液体介质密度，den2第二层液体介质密度;
//f源频率;
//x为待求参数，其他参数均作为已知值代入
void fadhrt2liq(double a,double b,double h,double eps,double x[],int m[],double H,double a1,double a2,double den1,double den2,double f)
{ 
	int n,js;
	double z,y,z1,y1,z0,y0;
	n=0;
	z=a;
	y=fadhrtf2liq(z,H,a1,a2,den1,den2,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=fadhrtf2liq(z,H,a1,a2,den1,den2,f);
		}
		
		else
		{  
			z1=z+h;y1=fadhrtf2liq(z1,H,a1,a2,den1,den2,f);      
			if(fabs(y1)<eps)
			{n=n+1;x[n-1]=z1;
			z=z1+h/2.0;y=fadhrtf2liq(z,H,a1,a2,den1,den2,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=fadhrtf2liq(z,H,a1,a2,den1,den2,f);
						js=1;
					}
					else
					{
						z0=(z1+z)/2.0;y0=fadhrtf2liq(z0,H,a1,a2,den1,den2,f);
						if(fabs(y0)<eps)
						{x[n]=z0;n=n+1;js=1;
						z=z0+h/2.0;
						y=fadhrtf2liq(z,H,a1,a2,den1,den2,f);
						}
						else if ((y*y0)<0.0)
						{z=z0;y=y0;} 
						else {z=z0;y=y0;}
					}
				}
			}
		}
	}
	m[1]=n;
	return;
}