discur3solid.h

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

#include"math.h"
//////////////////////////// fadhrtf3solidmo()三层固体介质频散函数算法实现////////////////////////////
//vr瑞雷波速度
//vs0,vs1,vs2分别为第一，二,三层的横波速度，vp0,vp1,vp2分别为第一,二,三层的纵波速度;
//h0,h1分别为第一，二层介质厚度;
//den0,den1,den2分别为第一，二，三层介质密度;
//f源频率;
double fadhrtf3solidmo(double vr,double vs0,double vs1,double vs2,double vp0,double vp1,double vp2,double h0,double h1,double den0,double den1,double den2,double f)
{	
    int i,j,i1;
	double cz=0.0,l[2],e[2],b[2],t[3],a[2],den[3],rp[3],g[3],vs[3],vp[3],d[2],r[2],s[2],h[2],rs[3],p[2],q[2];
	double	jz1[5][5],E1[5][1];		
	/////////////////////可调参数////////////////////////////////////
	vs[0]=vs0,vp[0]=vp0,h[0]=h0,den[0]=den0;
	vs[1]=vs1,vp[1]=vp1,h[1]=h1,den[1]=den1;
	vs[2]=vs2,vp[2]=vp2,den[2]=den2;



    t[0]=1.0-vr*vr/(2.0*vs[0]*vs[0]);
    t[1]=1.0-vr*vr/(2.0*vs[1]*vs[1]);
    t[2]=1.0-vr*vr/(2.0*vs[2]*vs[2]);
	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]);
	l[1]=den[2]*vs[2]*vs[2]/(den[1]*vs[1]*vs[1]);
	r[0]=vr*vr/(vp[0]*vp[0])-1.0;
	r[1]=vr*vr/(vp[1]*vp[1])-1.0;
    s[0]=vr*vr/(vs[0]*vs[0])-1.0;
    s[1]=vr*vr/(vs[1]*vs[1])-1.0;

	for(i=0;i<2;i++)
	{
		if(vr>vp[i])
		{
			rp[i]=sqrt(vr*vr/(vp[i]*vp[i])-1.0);
			p[i]=rp[i]*2.0*3.1415926*f*h[i]/vr;
			a[i]=cos(p[i]);
			e[i]=sin(p[i])/rp[i];

		}
		else
		{
			rp[i]=sqrt(1.0-vr*vr/(vp[i]*vp[i]));
			p[i]=rp[i]*2.0*3.1415926*f*h[i]/vr;
			a[i]=(exp(-1.0*p[i])+exp(p[i]))/2.0;
			e[i]=(exp(-1.0*p[i])-exp(p[i]))/(-2.0*rp[i]);			
		}
		///////////////////////////////////
		if(vr>vs[i])
		{
			rs[i]=sqrt(vr*vr/(vs[i]*vs[i])-1.0);	
			q[i]=rs[i]*2.0*3.1415926*f*h[i]/vr;
			b[i]=cos(q[i]);
			d[i]=sin(q[i])/rs[i];			          
		}
		else
		{
			rs[i]=sqrt(1.0-vr*vr/(vs[i]*vs[i]));
			q[i]=rs[i]*2.0*3.1415926*f*h[i]/vr;
			b[i]=(exp(-1.0*q[i])+exp(q[i]))/2.0;
			d[i]=(exp(-1.0*q[i])-exp(q[i]))/(-2.0*rs[i]);             
		}
		
		
	}	

	
    t[2]=1.0-vr*vr/(2.0*vs[2]*vs[2]);
	g[2]=1.0-t[2];
    rp[2]=sqrt(1.0-vr*vr/(vp[2]*vp[2]));
    rs[2]=sqrt(1.0-vr*vr/(vs[2]*vs[2]));
	
    	////////////////////////////////
	
	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 F2[5][5]={
		{(a[1]*b[1]*(1.0+t[1]*t[1])-e[1]*d[1]*(r[1]*s[1]+t[1]*t[1])-2.0*t[1])/l[1],2.0*(1.0-a[1]*b[1])*(1.0+t[1])+2.0*e[1]*d[1]*(r[1]+s[1]*t[1]),
		-1.0*(a[1]*d[1]+b[1]*e[1]*r[1])*g[1],(a[1]*d[1]*s[1]+b[1]*e[1])*g[1],(2.0*(1.0-a[1]*b[1])+e[1]*d[1]*(r[1]*s[1]+1.0))*l[1]},
		{((a[1]*b[1]-1.0)*(t[1]+t[1]*t[1])-e[1]*d[1]*(r[1]*s[1]+t[1]*t[1]*t[1]))/l[1],-4.0*a[1]*b[1]*t[1]+2.0*e[1]*d[1]*(r[1]*s[1]+t[1]*t[1])+(1.0+t[1])*(1.0+t[1]),
	    -1.0*(a[1]*d[1]*t[1]+b[1]*e[1]*r[1])*g[1],(a[1]*d[1]*s[1]+b[1]*e[1]*t[1])*g[1],((1.0-a[1]*b[1])*(1+t[1])+e[1]*d[1]*(r[1]*s[1]+t[1]))*l[1]},
		{(a[1]*d[1]*s[1]+b[1]*e[1]*t[1]*t[1])*g[1]/l[1],-2.0*(a[1]*d[1]*s[1]+b[1]*e[1]*t[1])*g[1],a[1]*b[1]*g[1]*g[1],e[1]*d[1]*s[1]*g[1]*g[1],-1.0*(a[1]*d[1]*s[1]+b[1]*e[1])*g[1]*l[1]},
		{-1.0*(a[1]*d[1]*t[1]*t[1]+b[1]*e[1]*r[1])*g[1]/l[1],2.0*(a[1]*d[1]*t[1]+b[1]*e[1]*r[1])*g[1],e[1]*d[1]*r[1]*g[1]*g[1],a[1]*b[1]*g[1]*g[1],(a[1]*d[1]+b[1]*e[1]*r[1])*g[1]*l[1]},
		{(2.0*(1.0-a[1]*b[1])*t[1]*t[1]+e[1]*d[1]*(r[1]*s[1]+pow(t[1],4)))/l[1],2.0*(a[1]*b[1]-1)*(t[1]+t[1]*t[1])-2.0*e[1]*d[1]*(r[1]*s[1]+t[1]*t[1]*t[1]),(a[1]*d[1]*t[1]*t[1]+b[1]*e[1]*r[1])*g[1],
    	-1.0*(a[1]*d[1]*s[1]+b[1]*e[1]*t[1]*t[1])*g[1],(a[1]*b[1]*(1.0+t[1]*t[1])-e[1]*d[1]*(r[1]*s[1]+t[1]*t[1])-2.0*t[1])*l[1]}
	};
	
	for(i=0;i<5;i++)
	{
		for(j=0;j<5;j++)
		{
			for(i1=0;i1<5;i1++) 								
				cz=F1[i][i1]*F2[i1][j]+cz;
			jz1[i][j]=cz;
			cz=0.0;
		}
		
	}
	
	
	
double	E3[5][1]={
		{1.0-rp[2]*rs[2]},
		{t[2]-rp[2]*rs[2]},
		{-1.0*rs[2]*g[2]},
		{rp[2]*g[2]},
		{-1.0*t[2]*t[2]+rp[2]*rs[2]}		
	};
	
	
	for(i=0;i<5;i++)
	{
		for(j=0;j<1;j++)
		{
			for(i1=0;i1<5;i1++) 								
				cz=jz1[i][i1]*E3[i1][j]+cz;
			E1[i][j]=cz;
			cz=0.0;
		}
		
	}
	return(E1[4][0]);
    
}
////////三层固体介质二分法计算///////////////
// fadhrt3solidmo()三层固体介质瑞雷波频散方程二分法求解算法,求解需用到函数fadhrtf3solidmo()
//a,b为根的求解区间，a<b，一般a大于三层介质中速度最小层介质的横波速度，b小于三层介质中速度最大层介质的横波速度;
//h为二分法求解步长;
//eps控制精度;
//x一维数组，长度为m[0]，返回在区间[a,b]中搜度到的瑞雷波速度的解;
//m一维数组，长度为2，一般取m[0]=10,取m[1]=0;
//vs0,vs1,vs2分别为第一，二,三层的横波速度，vp0,vp1,vp2分别为第一,二,三层的纵波速度;
//h0,h1分别为第一，二层介质厚度;
//den0,den1,den2分别为第一，二，三层介质密度;
//f源频率;
//x为待求参数，其他参数均作为已知值代入
void fadhrt3solidmo(double a,double b,double h,double eps,double x[],int m[],
					double vs0,double vs1,double vs2,double vp0,double vp1,double vp2,
					double h0,double h1,double den0,double den1,double den2,double f)
{ 
	int n,js;
	double z,y,z1,y1,z0,y0;
	n=0;
	z=a;
	y=fadhrtf3solidmo(z,vs0,vs1,vs2,vp0,vp1,vp2,h0,h1,den0,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=fadhrtf3solidmo(z,vs0,vs1,vs2,vp0,vp1,vp2,h0,h1,den0,den1,den2,f);
		}
		
		else
		{  
			z1=z+h;y1=fadhrtf3solidmo(z1,vs0,vs1,vs2,vp0,vp1,vp2,h0,h1,den0,den1,den2,f);      
			if(fabs(y1)<eps)
			{n=n+1;x[n-1]=z1;
			z=z1+h/2.0;y=fadhrtf3solidmo(z,vs0,vs1,vs2,vp0,vp1,vp2,h0,h1,den0,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=fadhrtf3solidmo(z,vs0,vs1,vs2,vp0,vp1,vp2,h0,h1,den0,den1,den2,f);
						js=1;
					}
					else
					{
						z0=(z1+z)/2.0;y0=fadhrtf3solidmo(z0,vs0,vs1,vs2,vp0,vp1,vp2,h0,h1,den0,den1,den2,f);
						if(fabs(y0)<eps)
						{x[n]=z0;n=n+1;js=1;
						z=z0+h/2.0;
						y=fadhrtf3solidmo(z,vs0,vs1,vs2,vp0,vp1,vp2,h0,h1,den0,den1,den2,f);
						}
						else if ((y*y0)<0.0)
						{z=z0;y=y0;} 
						else {z=z0;y=y0;}
					}
				}
			}
		}
	}
	m[1]=n;
	return;
}