VTI.zip_VTI_rugar_zip

clc clear all for n = 1:40; e = (n-1)* pi/180; e_grad(n) = n; vp = 3.30*1000; vs = 1.70*1000; D = 2.35*1000; vpp = 4.20*1000; vsp = 2.70*1000; Dp = 2.49*1000; f = asin((vs/vp)*sin(e)); ep = asin((vpp/vp)*sin(e)); fp = asin((vsp/vp)*sin(e)); A = [-sin(e) cos(f) sin(ep) cos(fp); cos(e) sin(f) cos(ep) sin(-fp); sin(2*e) -(vp/vs)*cos(2*f) (Dp/D)*(vp/vpp)*((vsp/vs)^2)*sin(2*ep) (Dp/D)*(vp/vsp)*((vsp/vs)^2)*cos(2*fp); -cos(2*f) -(vs/vp)*sin(2*f) (Dp/D)*(vpp/vp)*cos(2*fp) -(Dp/D)*(vsp/vp)*sin(2*fp)]; B = [sin(e); cos(e); sin(2*e); cos(2*f)]; X = inv(A)*B; A1(n) = (X(1)); del_vp = vpp-vp; vp_avg = (vp+vpp)/2; del_vs = vsp-vs; vs_avg = (vs +vsp)/2; del_D = Dp - D; D_avg = (D + Dp)/2; Z = D*vp; Zp = Dp*vpp; del_Z = Zp - Z; Z_avg = (Z+Zp)/2; sig = (0.5*((vp/vs)^2)-1)/(((vp/vs)^2)-1); sigp = (0.5*((vpp/vsp)^2)-1)/(((vpp/vsp)^2)-1); del_sig = sigp-sig; sig_avg = (sig + sigp)/2; G = D*((vs)^2); Gp = Dp*((vsp)^2); del_G = Gp-G; G_avg = (G+Gp)/2; delta =-0.120; epsilon = -0.133; f = asin((vs/vp)*sin(e)); ep = asin((vpp/vp)*sin(e)); fp = asin((vsp/vp)*sin(e)); A = [-sin(e) cos(f) sin(ep) cos(fp); cos(e) sin(f) cos(ep) sin(-fp); sin(2*e) -(vp/vs)*cos(2*f) (Dp/D)*(vp/vpp)*((vsp/vs)^2)*sin(2*ep) (Dp/D)*(vp/vsp)*((vsp/vs)^2)*cos(2*fp); -cos(2*f) -(vs/vp)*sin(2*f) (Dp/D)*(vpp/vp)*cos(2*fp) -(Dp/D)*(vsp/vp)*sin(2*fp)]; B = [sin(e); cos(e); sin(2*e); cos(2*f)]; X = inv(A)*B; A1(n) = (X(1)); %shuay approximation % sigma is poison ratio sigma = 0.319; sigmap = 0.148; sigma_avg = (sigma + sigmap)/2; % dZ/Z d_Z = 0.297; % dvp/vp d_vp = 0.240; % dD/D d_Dp = 0.058; A = 0.5*(d_Z); B = -2*(1 - sigma_avg/(1-sigma_avg))*A - 0.5*((1 - 3*sigma_avg)/(1 - sigma_avg))*d_vp + ((sigmap -sigma)/(1 - sigma_avg)^2); C = 0.5 * (d_vp); shuay(n) = A + B*(sin(e))^2 + C * (sin(e))^2* (tan (e))^2; % Rugar approximation % dG/G where G is shear impedance %d_G = 0.911 d_G = del_G/G_avg; vs_vpratio = 1.376; del_sigma = -0.172; Rugar(n) = 0.5 * (d_Z) + 0.5*(d_vp - (vs_vpratio * d_G))*(sin(e))^2 + 0.5 * d_vp *(sin(e))^2 * (tan(e))^2; %hilterman(1983) most used version of shuey equation hilterman(n) = (A * (cos(e))^2) + ( 2.25 *del_sigma * (sin(e))^2); %VTI approximation %approx_vti(n) = 0.5 * (d_Z) + 0.5*(d_vp - (vs_vpratio * d_G)+ delta)*(sin(e))^2 + 0.5 *( d_vp + epsilon) *(sin(e))^2 * (tan(e))^2; approx_vti(n) = 0.5 * (del_Z/Z_avg) + 0.5*((del_vp/vp_avg) - (((2*vs_avg)/vp_avg)^2 * (del_G/G_avg))+ delta )*(sin(e))^2 + 0.5 * ((del_vp/vp_avg)+epsilon) *(sin(e))^2 * (tan(e))^2; end for n = 1:40; e = (n-1)* pi/180; e_grad(n) = n; vp = 2.96*1000; vs = 1.38*1000; D = 2.43*1000; vpp = 3.49*1000; vsp = 2.29*1000; Dp = 2.14*1000; f = asin((vs/vp)*sin(e)); ep = asin((vpp/vp)*sin(e)); fp = asin((vsp/vp)*sin(e)); A = [-sin(e) cos(f) sin(ep) cos(fp); cos(e) sin(f) cos(ep) sin(-fp); sin(2*e) -(vp/vs)*cos(2*f) (Dp/D)*(vp/vpp)*((vsp/vs)^2)*sin(2*ep) (Dp/D)*(vp/vsp)*((vsp/vs)^2)*cos(2*fp); -cos(2*f) -(vs/vp)*sin(2*f) (Dp/D)*(vpp/vp)*cos(2*fp) -(Dp/D)*(vsp/vp)*sin(2*fp)]; B = [sin(e); cos(e); sin(2*e); cos(2*f)]; X = inv(A)*B; A1Z(n) = (X(1)); %shuay approximation % sigma is poison ratio sigma = 0.361; sigmap = 0.122; sigma_avg = (sigma + sigmap)/2; % dZ/Z d_Z = 0.038; % dvp/vp d_vp = 0.164; % dD/D d_Dp = -0.127; A = 0.5*(d_Z); B = -2*(1 - sigma_avg/(1-sigma_avg))*A - 0.5*((1 - 3*sigma_avg)/(1 - sigma_avg))*d_vp + ((sigmap -sigma)/(1 - sigma_avg)^2); C = 0.5 * (d_vp); shuayZ(n) = A + B*(sin(e))^2 + C * (sin(e))^2* (tan (e))^2; % Rugar approximation % dG/G where G is shear impedance d_G = 0.832; %(2vs/vp)^2 vs_vpratio = 1.295; del_sigma = -0.239; RugarZ(n) = 0.5 * (d_Z) + 0.5*(d_vp - (vs_vpratio * d_G))*(sin(e))^2 + 0.5 * d_vp *(sin(e))^2 * (tan(e))^2; %hilterman(1983) most used version of shuey equation hiltermanZ(n) = (A * (cos(e))^2) + ( 2.25 *del_sigma * (sin(e))^2); %VTI approximation approx_vtiZ(n) = 0.5 * (d_Z) + 0.5*(d_vp - (vs_vpratio * d_G)+ delta)*(sin(e))^2 + 0.5 *( d_vp + epsilon) *(sin(e))^2 * (tan(e))^2; end for n = 1:40; e = (n-1)* pi/180; e_grad(n) = n; vp = 2.73*1000; vs = 1.24*1000; D = 2.35*1000; vpp = 2.02*1000; vsp = 1.23*1000; Dp = 2.13*1000; f = asin((vs/vp)*sin(e)); ep = asin((vpp/vp)*sin(e)); fp = asin((vsp/vp)*sin(e)); A = [-sin(e) cos(f) sin(ep) cos(fp); cos(e) sin(f) cos(ep) sin(-fp); sin(2*e) -(vp/vs)*cos(2*f) (Dp/D)*(vp/vpp)*((vsp/vs)^2)*sin(2*ep) (Dp/D)*(vp/vsp)*((vsp/vs)^2)*cos(2*fp); -cos(2*f) -(vs/vp)*sin(2*f) (Dp/D)*(vpp/vp)*cos(2*fp) -(Dp/D)*(vsp/vp)*sin(2*fp)]; B = [sin(e); cos(e); sin(2*e); cos(2*f)]; X = inv(A)*B; A1N(n) = (X(1)); %shuay approximation % sigma is poison ratio sigma = 0.370; sigmap = 0.205; sigma_avg = (sigma + sigmap)/2; % dZ/Z d_Z = -0.394; % dvp/vp d_vp = -0.298; % dD/D d_Dp = -0.098; A = 0.5*(d_Z); B = -2*(1 - sigma_avg/(1-sigma_avg))*A - 0.5*((1 - 3*sigma_avg)/(1 - sigma_avg))*d_vp + ((sigmap -sigma)/(1 - sigma_avg)^2); C = 0.5 * (d_vp); shuayN(n) = A + B*(sin(e))^2 + C * (sin(e))^2* (tan (e))^2; % Rugar approximation % dG/G where G is shear impedance d_G = -0.114; %(2vs/vp)^2 vs_vpratio = 1.081; del_sigma = -0.164; RugarN(n) = 0.5 * (d_Z) + 0.5*(d_vp - (vs_vpratio * d_G))*(sin(e))^2 + 0.5 * d_vp *(sin(e))^2 * (tan(e))^2; %hilterman(1983) most used version of shuey equation hiltermanN(n) = (A * (cos(e))^2) + ( 2.25 *del_sigma * (sin(e))^2); %VTI approximation approx_vtiN(n) = 0.5 * (d_Z) + 0.5*(d_vp - (vs_vpratio * d_G)+ delta)*(sin(e))^2 + 0.5 *( d_vp + epsilon) *(sin(e))^2 * (tan(e))^2; end % figures figure(1) plot(e_grad,A1,'b','linewidth',2); hold on plot (e_grad,approx_vti, 'g','linewidth',2); hold on legend({'exact','approx-VTI'},'fontsize',20) %plot (e_grad, shuay,'k','linewidth',2); %hold on %plot (e_grad,Rugar,'g','linewidth',2); %hold on plot(e_grad,A1Z,'b','linewidth',2); %set(graph4,'linewidth',2) hold on %graph5=plot(e_grad,shuayZ,'k'); %set(graph5,'linewidth',2) %hold on %graph6=plot(e_grad,RugarZ,'g'); %set(graph6,'linewidth',2) %hold on plot(e_grad,A1N,'b','linewidth',2); hold on %set(graph7,'linewidth',2) %hold on %graph8=plot(e_grad,shuayN,'k'); %set(graph8,'linewidth',2) %hold on %graph9=plot(e_grad,RugarN,'g'); %xlabel('incidence angle (Deg)','fontsize',20); %ylabel('P-Wave RC','fontsize',20); %legend ({'exact','shuay','Rugar'},'fontsize',20) %figure(2) %plot(e_grad,hilterman,'m','linewidth',2) %hold on %plot (e_grad, shuay,'k','linewidth',2); %hold on %plot(e_grad,hiltermanZ,'m','linewidth',2) %hold on %graph2 =plot (e_grad, shuayZ,'k','linewidth',2); %hold on %plot(e_grad,hiltermanN,'m','linewidth',2) %hold on %graph2 =plot (e_grad, shuayN,'k','linewidth',2 ); %legend ('shuay','hilterman'); hold on hold on plot (e_grad,approx_vtiZ, 'g','linewidth',2); hold on plot (e_grad,approx_vtiN, 'g','linewidth',2); xlabel('incidence angle (Deg)','fontsize',20); ylabel('P-Wave RC','fontsize',20); title('VTI Rugar approximation'); grid on