functionsps.h

实现共轭梯度法的实例

#include "math.h"
float Rps_Aki(float IncAngle,float Vp1,float Vp2,float Vs1,float Vs2,float Rho1,float Rho2,float& d,float& e)
//注：这里的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 fai1=asin(Vs1*sin(sita1)/Vp1);
	float fai2=asin(Vs2*sin(sita1)/Vp1);
	float fai=0.5*(fai1+fai2);

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

	D=(0.5*alpha*tan(fai)/beta)*
	(4*pow(sin(fai),2)-4*(beta/alpha)*cos(sita)*cos(fai));
	
	
	E=-(2*alpha*tan(fai)/beta)*(1+2*pow(sin(fai),2)-2*(beta/alpha)*cos(sita)*cos(fai));
	

	
	item1=D*d_J/J;
	item2=E*d_Rho/Rho;
	
	r_ps=item1+item2;
	d=D;
	e=E;
	return r_ps;
}