discur2solid.h

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

////////////////////////////两层固体介质无量纲快速矢量算法///////////////////////////////
#include"math.h"
////////////////////////////fadhrtfdoublesoliddiscur()频散函数算法实现////////////////////////////
//vr瑞雷波速度
//vs0,vs1分别为第一，二层介质的横波速度
//vp0,vp1分别为第一，二层介质的纵波速度
//h0第一层介质厚度
//den0,den1分别为第一，二层介质密度
//f源频率
double fadhrtfdoublesoliddiscur(double vr,double vs0,double vs1,double vp0,double vp1,double h0, double den0,double den1,double f )
{
	
    int i,j,i1;
	double cz=0.0,l[1],e[1],b[1],t[2],a[1],den[2],rp[2],g[2],vs[2],vp[2],d[1],r[1],s[1],h[1],rs[2],p[1],q[1];
	double	E1[5][1];		
	/////////////////////可调参数////////////////////////////////////
	vs[0]=vs0,vp[0]=vp0,h[0]=h0,den[0]=den0;
	vs[1]=vs1,vp[1]=vp1,den[1]=den1;	
	rp[0]=sqrt(1-vr*vr/pow(vp[0],2));
	p[0]=rp[0]*2*3.1415926*f*h[0]/vr;
	a[0]=(exp(-1.0*p[0])+exp(p[0]))/2.0;
	e[0]=(exp(-1.0*p[0])-exp(p[0]))/(-2.0*rp[0]);
	rs[0]=sqrt(vr*vr/(vs[0]*vs[0])-1.0);
	q[0]=rs[0]*2.0*3.1415926*f*h[0]/vr;
	b[0]=cos(q[0]);
	d[0]=sin(q[0])/rs[0];	
	rp[1]=sqrt(1.0-vr*vr/(vp[1]*vp[1]));	
	rs[1]=sqrt(1.0-vr*vr/(vs[1]*vs[1]));	
    t[0]=1.0-vr*vr/(2.0*vs[0]*vs[0]);
    t[1]=1.0-vr*vr/(2.0*vs[1]*vs[1]);
    r[0]=vr*vr/pow(vp[0],2)-1.0;
	s[0]=vr*vr/(vs[0]*vs[0])-1.0;	
    g[0]=1.0-t[0];
    g[1]=1.0-t[1];
	l[0]=den[1]*vs[1]*vs[1]/(den[0]*vs[0]*vs[0]);	
	double F1[5][5]={
		{(a[0]*b[0]*(1.0+t[0]*t[0])-e[0]*d[0]*(r[0]*s[0]+t[0]*t[0])-2.0*t[0])/l[0],2.0*(1.0-a[0]*b[0])*(1.0+t[0])+2.0*e[0]*d[0]*(r[0]+s[0]*t[0]),
			-1.0*(a[0]*d[0]+b[0]*e[0]*r[0])*g[0],(a[0]*d[0]*s[0]+b[0]*e[0])*g[0],(2.0*(1.0-a[0]*b[0])+e[0]*d[0]*(r[0]*s[0]+1.0))*l[0]},
		{((a[0]*b[0]-1.0)*(t[0]+t[0]*t[0])-e[0]*d[0]*(r[0]*s[0]+t[0]*t[0]*t[0]))/l[0],-4.0*a[0]*b[0]*t[0]+2.0*e[0]*d[0]*(r[0]*s[0]+t[0]*t[0])+(1.0+t[0])*(1.0+t[0]),
		-1.0*(a[0]*d[0]*t[0]+b[0]*e[0]*r[0])*g[0],(a[0]*d[0]*s[0]+b[0]*e[0]*t[0])*g[0],((1.0-a[0]*b[0])*(1+t[0])+e[0]*d[0]*(r[0]*s[0]+t[0]))*l[0]},
		{(a[0]*d[0]*s[0]+b[0]*e[0]*t[0]*t[0])*g[0]/l[0],-2.0*(a[0]*d[0]*s[0]+b[0]*e[0]*t[0])*g[0],a[0]*b[0]*g[0]*g[0],e[0]*d[0]*s[0]*g[0]*g[0],-1.0*(a[0]*d[0]*s[0]+b[0]*e[0])*g[0]*l[0]},
		{-1.0*(a[0]*d[0]*t[0]*t[0]+b[0]*e[0]*r[0])*g[0]/l[0],2.0*(a[0]*d[0]*t[0]+b[0]*e[0]*r[0])*g[0],e[0]*d[0]*r[0]*g[0]*g[0],a[0]*b[0]*g[0]*g[0],(a[0]*d[0]+b[0]*e[0]*r[0])*g[0]*l[0]},
		{(2.0*(1.0-a[0]*b[0])*t[0]*t[0]+e[0]*d[0]*(r[0]*s[0]+pow(t[0],4)))/l[0],2.0*(a[0]*b[0]-1)*(t[0]+t[0]*t[0])-2.0*e[0]*d[0]*(r[0]*s[0]+t[0]*t[0]*t[0]),(a[0]*d[0]*t[0]*t[0]+b[0]*e[0]*r[0])*g[0],
		-1.0*(a[0]*d[0]*s[0]+b[0]*e[0]*t[0]*t[0])*g[0],(a[0]*b[0]*(1.0+t[0]*t[0])-e[0]*d[0]*(r[0]*s[0]+t[0]*t[0])-2.0*t[0])*l[0]}
	};	
	double	E3[5][1]={
		{1.0-rp[1]*rs[1]},
		{t[1]-rp[1]*rs[1]},
		{-1.0*rs[1]*g[1]},
		{rp[1]*g[1]},
		{-1.0*t[1]*t[1]+rp[1]*rs[1]}		
	};
	
	for(i=0;i<5;i++)
	{
		for(j=0;j<1;j++)
		{
			for(i1=0;i1<5;i1++) 								
				cz=F1[i][i1]*E3[i1][j]+cz;
			E1[i][j]=cz;
			cz=0.0;
		}		
	}
	return(E1[4][0]);  
}
//fadhrt2solid()两层固体介质频散方程二分法求解算法,求解需用到函数fadhrtfdoublesoliddiscur()
//a,b为根的求解区间，a<b，a大于速度较小固体层横波速度，b小于速度较大固体层横波速度;
//h为二分法求解步长;
//eps控制精度;
//x一维数组，长度为m[0]，返回在区间[a,b]中搜度到的瑞雷波速度的解;
//m一维数组，长度为2，一般取m[0]=10,取m[1]=0;
//vs0,vs1分别为第一，二层介质的横波速度
//vp0,vp1分别为第一，二层介质的纵波速度
//h0第一层介质厚度
//den0,den1分别为第一，二层介质密度
//f源频率
void fadhrt2solid(double a1,double b,double h,double eps,double x[],int m[],double vs0,double vs1,double vp0,double vp1,
double h0, double den0,double den1,double f)
		{ 
			int n,js;
			double z,y,z1,y1,z0,y0;
			n=0;
			z=a1;
			y=fadhrtfdoublesoliddiscur(z,vs0,vs1,vp0,vp1,h0,den0,den1,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=fadhrtfdoublesoliddiscur(z,vs0,vs1,vp0,vp1,h0,den0,den1,f);
				}
				
				else
				{  
					z1=z+h;y1=fadhrtfdoublesoliddiscur(z1,vs0,vs1,vp0,vp1,h0,den0,den1,f);      
					if(fabs(y1)<eps)
					{n=n+1;x[n-1]=z1;
					z=z1+h/2.0;y=fadhrtfdoublesoliddiscur(z,vs0,vs1,vp0,vp1,h0,den0,den1,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=fadhrtfdoublesoliddiscur(z,vs0,vs1,vp0,vp1,h0,den0,den1,f);
								js=1;
							}
							else
							{
								z0=(z1+z)/2.0;y0=fadhrtfdoublesoliddiscur(z0,vs0,vs1,vp0,vp1,h0,den0,den1,f);
								if(fabs(y0)<eps)
								{x[n]=z0;n=n+1;js=1;
								z=z0+h/2.0;
								y=fadhrtfdoublesoliddiscur(z,vs0,vs1,vp0,vp1,h0,den0,den1,f);
								}
								else if ((y*y0)<0.0)
								{z=z0;y=y0;} 
								else {z=z0;y=y0;}
							}
						}
					}
				}
			}
			m[1]=n;
			return;
		}