functions.h

实现共轭梯度法的实例

#include "math.h"
float Rpp_Aki(float IncAngle,float Vp1,float Vp2,float Vs1,float Vs2,float Rho1,float Rho2)
//注：这里的sita为sita1即入射角
{
	float alpha=(Vp1+Vp2)/2;
	float beta= (Vs1+Vs2)/2;
	float Rho = (Rho1+Rho2)/2;
	float I=(Vp1*Rho1+Vp2*Rho2)/2;
	float J=(Vs1*Rho1+Vs2*Rho2)/2;

	float d_alpha=Vp2-Vp1;
	float d_beta =Vs2-Vs1;
	float d_Rho  =Rho2-Rho1;
	float d_I=Vp2*Rho2-Vp1*Rho1;
	float d_J=Vs2*Rho2-Vs1*Rho1;
	float sita1=0,sita2=0,sita=0;
	sita1=IncAngle;
	sita2=asin(Vp2*sin(sita1)/Vp1);
	sita=0.5*(sita1+sita2);

	
	float r_pp=0;
	float A=0,B=0,C=0;   //分别为纵波阻抗、横波阻抗以及密度差异前的系数
	float item1=0,item2=0,item3=0;//分别为三个参数对应的三项

	A=0.5*(1+pow(tan(sita),2));
	B=-4*pow(beta*sin(sita)/alpha,2);
	C=2*pow(beta*sin(sita)/alpha,2)-0.5*pow(tan(sita),2);

	item1=A*d_I/I;;
	item2=B*d_J/J;
	item3=C*d_Rho/Rho;
	
	r_pp=item1+item2+item3;
	return r_pp;
}

float Rpp_Bortfeld(float IncAngle,float Vp1,float Vp2,float Vs1,float Vs2,float Rho1,float Rho2)
{
	float item1=0,item2=0;
	float sita1=0,sita2=0;
	sita1=IncAngle;
	sita2=asin(Vp2*sin(sita1)/Vp1);
	float r_pp=0;
	item1=0.5*log(Vp2*Rho2*cos(sita1)/(Vp1*Rho1*cos(sita2)));
	item2=pow(sin(sita1)/Vp1,2)*(pow(Vs1,2)-pow(Vs2,2))*(2+ log(Rho2/Rho1)/(log(Vp2/Vp1)-log(Vp2*Vs2/(Vp1*Vs1)) ));
	
	r_pp=item1+item2;
	return r_pp;
}


float Rpp_Zheng(float IncAngle,float Vp1,float Vp2,float Vs1,float Vs2,float Rho1,float Rho2)
{
	float R0=0,R2=0,R4=0;

	float alpha=(Vp1+Vp2)/2;
	float beta= (Vs1+Vs2)/2;
	float Rho = (Rho1+Rho2)/2;

	float d_alpha=Vp2-Vp1;
	float d_beta =Vs2-Vs1;
	float d_Rho  =Rho2-Rho1;
	float r_pp;

	float sita1=0,sita2=0,sita=0;
	sita1=IncAngle;
	sita2=asin(Vp2*sin(sita1)/Vp1);
	sita=0.5*(sita1+sita2);

	R0=0.5*(d_alpha/alpha+d_Rho/Rho);
	R2=0.5*d_alpha/alpha
	   -2*pow(beta/alpha,2)*d_Rho/Rho
	   -4*pow(beta/alpha,2)*d_beta/beta;
	R4=0.5*d_alpha/alpha;
	r_pp=R0+R2*pow(sin(sita),2)+R4*pow(sin(sita),4);
	return r_pp;
}


float Rpp_Cerveny(float IncAngle,float Vp1,float Vp2,float Vs1,float Vs2,float Rho1,float Rho2)
{
	
	float Nr=0;
	float D=0;
	float r_pp;

	float p=sin(IncAngle)/Vp1;
	float p2=p*p;

	float Vp12=Vp1*Vp1;
	float Vs12=Vs1*Vs1;
	float Vp22=Vp2*Vp2;
	float Vs22=Vs2*Vs2;

	float q=2*(Rho2*Vs22-Rho1*Vs12);
	float q2=q*q;
	

	float X=Rho2-q*p2;
	float Y=Rho1+q*p2;
	float Z=Rho2-Rho1-q*p2;

	float P1=pow(1-Vp12*p2,0.5);
	float P2=pow(1-Vs12*p2,0.5);
	float P3=pow(1-Vp22*p2,0.5);
	float P4=pow(1-Vs22*p2,0.5);
	


	float X2=X*X;
	float Y2=Y*Y;
	float Z2=Z*Z;

	Nr=q2*p2*P1*P2*P3*P4
	   +Rho1*Rho2*(Vs1*Vp2*P1*P4-Vp1*Vs2*P2*P3)
	   -Vp1*Vs1*P3*P4*Y2
	   +Vp2*Vs2*P1*P2*X2
	   -Vp1*Vp2*Vs1*Vs2*p2*Z2;

	D=q2*p2*P1*P2*P3*P4
	   +Rho1*Rho2*(Vs1*Vp2*P1*P4+Vp1*Vs2*P2*P3)
	   +Vp1*Vs1*P3*P4*Y2
	   +Vp2*Vs2*P1*P2*X2
	   +Vp1*Vp2*Vs1*Vs2*p2*Z2;

    r_pp=Nr/D;

	return r_pp;
}

float Rpp_Exact(float IncAngle,float Vp1,float Vp2,float Vs1,float Vs2,float Rho1,float Rho2)
{
	float r_pp;
	float i1=IncAngle;
	float p=sin(i1)/Vp1;
	float i2=asin(Vp2*p);
	float j1=asin(Vs1*p);
	float j2=asin(Vs2*p);

	float p2=p*p;
	float beta12=Vs1*Vs1;
	float beta22=Vs2*Vs2;

	float a=Rho2*(1-2*beta22*p2)-Rho1*(1-2*beta12*p2);
	float b=Rho2*(1-2*beta22*p2)+2*Rho1*beta12*p2;
	float c=Rho1*(1-2*beta12*p2)+2*Rho2*beta22*p2;
	float d=2*(Rho2*beta22-Rho1*beta12);

	float E=  b*cos(i1)/Vp1+c*cos(i2)/Vp2;
	float F=  b*cos(j1)/Vs1+c*cos(j2)/Vs2;
	float G=a-d*cos(i1)*cos(j2)/(Vp1*Vs2);
	float H=a-d*cos(i2)*cos(j1)/(Vp2*Vs1);
	float D=E*F+G*H*p2;
	
	r_pp=((b*cos(i1)/Vp1-c*cos(i2)/Vp2)*F-(a+d*cos(i1)*cos(j2)/(Vp1*Vp2))*H*p2)/D;
	return r_pp;
}