discur4solid.h

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

//////////////////////////四层固体介质频散函数计算///////////////////////////////////
#include"math.h" 
//////////////////////////// fadhrtf4solidmodel()四层固体层状介质频散函数算法实现////////////////////////////
//vr瑞雷波速度；
//vs0,vs1,vs2,vs3分别为第一，二,三,四层的横波速度，vp0,vp1,vp2,vp3分别为第一,二,三,四层的纵波速度;
//h0,h1,h2分别为第一，二,三层介质厚度;
//den0,den1,den2，den3分别为第一，二，三,四层介质密度;
//f源频率;
double fadhrtf4solidmodel(double vr,double vs0,double vs1,double vs2,double vs3,double vp0,double vp1,double vp2,
			  double vp3,double h0,double h1,double h2,double den0,double den1,double den2,double den3,double f)
{
	
    int i,j,i1;
	double cz=0.0,l[3],e[3],b[3],t[4],a[3],den[4],rp[4],g[4],vs[4],vp[4],d[3],r[3],s[3],h[3],rs[4],p[3],q[3];
	double	jz1[5][5],jz2[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,h[2]=h2,den[2]=den2;
    vs[3]=vs3,vp[3]=vp3,den[3]=den3;
	
	
	 rp[3]=sqrt(1.0-vr*vr/(vp[3]*vp[3]));
     rs[3]=sqrt(1.0-vr*vr/(vs[3]*vs[3]));
     t[3]=1.0-vr*vr/(2.0*vs[3]*vs[3]);
	 g[3]=1.0-t[3]; 

	for(i=0;i<3;i++)
	{
		t[i]=1.0-vr*vr/(2.0*vs[i]*vs[i]);
		g[i]=1.0-t[i];
		l[i]=den[i+1]*vs[i+1]*vs[i+1]/(den[i]*vs[i]*vs[i]);
		r[i]=vr*vr/(vp[i]*vp[i])-1.0;
        s[i]=vr*vr/(vs[i]*vs[i])-1.0;

		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[3]=1.0-vr*vr/(2.0*vs[3]*vs[3]);
	g[3]=1.0-t[3];
    rp[3]=sqrt(1.0-vr*vr/(vp[3]*vp[3]));
    rs[3]=sqrt(1.0-vr*vr/(vs[3]*vs[3]));
	
    	////////////////////////////////
	
	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]}
	};

		double F3[5][5]={
		{(a[2]*b[2]*(1.0+t[2]*t[2])-e[2]*d[2]*(r[2]*s[2]+t[2]*t[2])-2.0*t[2])/l[2],2.0*(1.0-a[2]*b[2])*(1.0+t[2])+2.0*e[2]*d[2]*(r[2]+s[2]*t[2]),
		-1.0*(a[2]*d[2]+b[2]*e[2]*r[2])*g[2],(a[2]*d[2]*s[2]+b[2]*e[2])*g[2],(2.0*(1.0-a[2]*b[2])+e[2]*d[2]*(r[2]*s[2]+1.0))*l[2]},
		{((a[2]*b[2]-1.0)*(t[2]+t[2]*t[2])-e[2]*d[2]*(r[2]*s[2]+t[2]*t[2]*t[2]))/l[2],-4.0*a[2]*b[2]*t[2]+2.0*e[2]*d[2]*(r[2]*s[2]+t[2]*t[2])+(1.0+t[2])*(1.0+t[2]),
	    -1.0*(a[2]*d[2]*t[2]+b[2]*e[2]*r[2])*g[2],(a[2]*d[2]*s[2]+b[2]*e[2]*t[2])*g[2],((1.0-a[2]*b[2])*(1+t[2])+e[2]*d[2]*(r[2]*s[2]+t[2]))*l[2]},
		{(a[2]*d[2]*s[2]+b[2]*e[2]*t[2]*t[2])*g[2]/l[2],-2.0*(a[2]*d[2]*s[2]+b[2]*e[2]*t[2])*g[2],a[2]*b[2]*g[2]*g[2],e[2]*d[2]*s[2]*g[2]*g[2],-1.0*(a[2]*d[2]*s[2]+b[2]*e[2])*g[2]*l[2]},
		{-1.0*(a[2]*d[2]*t[2]*t[2]+b[2]*e[2]*r[2])*g[2]/l[2],2.0*(a[2]*d[2]*t[2]+b[2]*e[2]*r[2])*g[2],e[2]*d[2]*r[2]*g[2]*g[2],a[2]*b[2]*g[2]*g[2],(a[2]*d[2]+b[2]*e[2]*r[2])*g[2]*l[2]},
		{(2.0*(1.0-a[2]*b[2])*t[2]*t[2]+e[2]*d[2]*(r[2]*s[2]+pow(t[2],4)))/l[2],2.0*(a[2]*b[2]-1)*(t[2]+t[2]*t[2])-2.0*e[2]*d[2]*(r[2]*s[2]+t[2]*t[2]*t[2]),(a[2]*d[2]*t[2]*t[2]+b[2]*e[2]*r[2])*g[2],
    	-1.0*(a[2]*d[2]*s[2]+b[2]*e[2]*t[2]*t[2])*g[2],(a[2]*b[2]*(1.0+t[2]*t[2])-e[2]*d[2]*(r[2]*s[2]+t[2]*t[2])-2.0*t[2])*l[2]}
	};
	
	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;
		}
		
	}
	
    	for(i=0;i<5;i++)
	{
		for(j=0;j<5;j++)
		{
			for(i1=0;i1<5;i1++) 								
				cz=jz1[i][i1]*F3[i1][j]+cz;
			jz2[i][j]=cz;
			cz=0.0;
		}
		
	}	
	
double	E4[5][1]={
		{1.0-rp[3]*rs[3]},
		{t[2]-rp[3]*rs[3]},
		{-1.0*rs[3]*g[3]},
		{rp[3]*g[3]},
		{-1.0*t[3]*t[3]+rp[3]*rs[3]}		
	};
	
	
	for(i=0;i<5;i++)
	{
		for(j=0;j<1;j++)
		{
			for(i1=0;i1<5;i1++) 								
				cz=jz2[i][i1]*E4[i1][j]+cz;
			E1[i][j]=cz;
			cz=0.0;
		}
		
	}
	return(E1[4][0]);
    
}

//////////////////////////fadhrt4solidmodel()二分法求四层固体介质频散曲线///////////////////////////////////
// fadhrt4solidmodel()四层固体介质频散方程二分法求解算法,求解需用到函数fadhrtf4solidmodel()
//a,b为根的求解区间，a<b，一般a大于四层介质中速度最小层介质的横波速度，b小于四层介质中速度最大层介质的横波速度;
//h为二分法求解步长;
//eps控制精度;
//x一维数组，长度为m[0]，返回在区间[a,b]中搜度到的瑞雷波速度的解;
//m一维数组，长度为2，一般取m[0]=10,取m[1]=0;
//vs0,vs1,vs2，vs3分别为第一，二,三,四层的横波速度，vp0,vp1,vp2,vp3分别为第一,二,三，四层的纵波速度;
//h0,h1,h2分别为第一，二,三层介质厚度;
//den0,den1,den2,den3分别为第一，二，三,四层介质密度;
//f源频率;
//x为待求参数，其他参数均作为已知值代入
	void fadhrt4solidmodel(double a1,double b,double h,double eps,double x[],
	int m[],double vs0,double vs1,double vs2,
	double vs3,double vp0,double vp1,double vp2,
			  double vp3,double h0,double h1,double h2,double den0,
			  double den1,double den2,double den3,double f)
		{ 
			int n,js;
			double z,y,z1,y1,z0,y0;
			n=0;
			z=a1;
			y=fadhrtf4solidmodel(z,vs0,vs1,vs2,vs3,vp0,vp1,vp2,vp3,h0,h1,h2,den0,den1,den2,den3,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=fadhrtf4solidmodel(z,vs0,vs1,vs2,vs3,vp0,vp1,vp2,vp3,h0,h1,h2,den0,den1,den2,den3,f);
				}
				
				else
				{  
					z1=z+h;y1=fadhrtf4solidmodel(z1,vs0,vs1,vs2,vs3,vp0,vp1,vp2,vp3,h0,h1,h2,den0,den1,den2,den3,f);      
					if(fabs(y1)<eps)
					{n=n+1;x[n-1]=z1;
					z=z1+h/2.0;y=fadhrtf4solidmodel(z,vs0,vs1,vs2,vs3,vp0,vp1,vp2,vp3,h0,h1,h2,den0,den1,den2,den3,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=fadhrtf4solidmodel(z,vs0,vs1,vs2,vs3,vp0,vp1,vp2,vp3,h0,h1,h2,den0,den1,den2,den3,f);
								js=1;
							}
							else
							{
								z0=(z1+z)/2.0;y0=fadhrtf4solidmodel(z0,vs0,vs1,vs2,vs3,vp0,vp1,vp2,vp3,h0,h1,h2,den0,den1,den2,den3,f);
								if(fabs(y0)<eps)
								{x[n]=z0;n=n+1;js=1;
								z=z0+h/2.0;
								y=fadhrtf4solidmodel(z,vs0,vs1,vs2,vs3,vp0,vp1,vp2,vp3,h0,h1,h2,den0,den1,den2,den3,f);
								}
								else if ((y*y0)<0.0)
								{z=z0;y=y0;} 
								else {z=z0;y=y0;}
							}
						}
					}
				}
			}
			m[1]=n;
			return;
		}