liqsolid.h

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

///////////////////计算二层介质（上层为液体，下层为固体）面波频散函数//////////////////////
#include"math.h"
////////////////////////////fadhrtf2liqsolid()频散函数算法实现////////////////////////////
//vr瑞雷波速度
//a1液体介质纵波速度
//a2固体介质纵波速度
//b2固体介质横波速度
//den1液体介质密度,den2固体介质密度
//f源频率
//H液体介质厚度
double fadhrtf2liqsolid(double vr,double a1,double a2,double b2,double den1,double den2,double f,double H)
{
     double zuo,you,out;
    ///////////////////////////////////////////////kkk/////////////////////////
    zuo=sqrt(1.0-vr*vr/pow(a2,2))*tan((2.0*3.1415926*f/vr)*H*sqrt(vr*vr/pow(a1,2)-1.0));
	you=(den2*pow(b2,4)/(den1*pow(vr,4)))*sqrt(vr*vr/(a1*a1)-1.0)*(4.0*sqrt(1.0-vr*vr/(a2*a2))*sqrt(1.0-vr*vr/(b2*b2))-pow(2.0-vr*vr/(b2*b2),2));
	out=zuo-you;		
	return(out);	
}
////两层，上层液体，下层固体二分法程序////////
//fadhrt2liqsolid()四层固体介质频散方程二分法求解算法,求解需用到函数fadhrtf2liqsolid()
//a,b为根的求解区间，a<b，一般a大于液体介质横波速度，b小于固体介质横波速度;
//h为二分法求解步长;
//eps控制精度;
//x一维数组，长度为m[0]，返回在区间[a,b]中搜度到的瑞雷波速度的解;
//m一维数组，长度为2，一般取m[0]=10,取m[1]=0;
//a1液体介质纵波速度
//a2固体介质纵波速度
//b2固体介质横波速度
//den1液体介质密度,den2固体介质密度
//H液体层介质厚度
//f源频率
//x为待求参数，其他参数均作为已知值代入
void fadhrt2liqsolid(double a,double b,double h,double eps,double x[],int m[],double a1,double a2,
			double b2,double den1,double den2,double H,double f)
{ 
	int n,js;
	double z,y,z1,y1,z0,y0;
	n=0;
	z=a;
	y=fadhrtf2liqsolid(z,a1,a2,b2,den1,den2,H,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=fadhrtf2liqsolid(z,a1,a2,b2,den1,den2,H,f);
		}
		
		else
		{  
			z1=z+h;y1=fadhrtf2liqsolid(z1,a1,a2,b2,den1,den2,H,f);      
			if(fabs(y1)<eps)
			{n=n+1;x[n-1]=z1;
			z=z1+h/2.0;y=fadhrtf2liqsolid(z,a1,a2,b2,den1,den2,H,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=fadhrtf2liqsolid(z,a1,a2,b2,den1,den2,H,f);
						js=1;
					}
					else
					{
						z0=(z1+z)/2.0;y0=fadhrtf2liqsolid(z0,a1,a2,b2,den1,den2,H,f);
						if(fabs(y0)<eps)
						{x[n]=z0;n=n+1;js=1;
						z=z0+h/2.0;
						y=fadhrtf2liqsolid(z,a1,a2,b2,den1,den2,H,f);
						}
						else if ((y*y0)<0.0)
						{z=z0;y=y0;} 
						else {z=z0;y=y0;}
					}
				}
			}
		}
	}
	m[1]=n;
	return;
}