www.pudn.com > avo_shuey.rar > avo_shuey.m, change:2010-05-12,size:1005b


function [RT,sinp]=avo_shuey(density1,density2,Vp1,Vp2,sigma1,sigma2)
%sigma:泊松比
p=0:pi/180:pi/2; %入射角
sigma2=[0.0 0.1 0.2 0.3 0.4];
Vp=(Vp1+Vp2)/2;
dVp=Vp2-Vp1;
density=(density1+density2)/2;
ddensity=density2-density1;
R0=0.5*(dVp/Vp+ddensity/density);
B=(dVp/Vp)/(dVp/Vp+ddensity/density);
for i=1:length(sigma2)
    sigma(i)=(sigma1+sigma2(i))/2;
    dsigma(i)=sigma2(i)-sigma1;
%     Vs=Vp*sqrt(0.5*(1-2*sigma(i))/(1-sigma));
    A0(i)=B-2*(1+B)*(1-2*sigma(i))/(1-sigma(i));
    A(i)=A0(i)+1/(1-sigma(i))^2*dsigma(i)/R0;
    for j=1:length(p)
        %     RT(i)=R0+(A0*R0+dsigma/(1-sigma)^2)*sin(p(i))^2+0.5*dVp/Vp*(tan(p(i))^2-sin(p(i))^2);%Shuey完全形式
        RR(j)=R0+(A0(i)*R0+dsigma(i)/(1-sigma(i))^2)*sin(p(j))^2;%Shuey略去高次项
%         RR(j)=R0*(1+A(i)*p(j)^2);
        sinp(j)=sin(p(j))^2;
    end
    RT(i,:)=RR;
end

% i=0:30;
% figure(1)  %纵波反射系数图
% plot(i,RT(1,:),'b');
% xlabel('入射角(°)');ylabel('纵波反射系数');title('反射P波');
% axis([0 30 -1 2]);
% grid on;