%REFLPSL calcule coef. reflexion-transmission P-SV interface liquide-solide
% voir : http://jmmeost.free.fr/Sismique/osis.html#PSI
% fig1  = http://jmmeost.free.fr/Sismique/reflpsl.gif

clear all, clf

% cas sediment mou
VP1 = 1500 ; VP2 = 1700; VS2 = VP2/2; 
ro1 = 1; ro2 = 1.6;

deg = 180/pi;
th1c = asin(VP1/VP2); % angle critique

th1 = .001:.002:pi/2; % angle incidence dans eau
p = sin(th1)/VP1;

a = -ro2*sqrt(1-(p*VP1).^2)/VP1.*((1-2*(p*VS2).^2).^2+4*VS2^4*p.^2.*sqrt(1-(p*VP2).^2)/VP2.*sqrt(1-(p*VS2).^2)/VS2);
b = -ro1*sqrt(1-(p*VP2).^2)/VP2;

DL = a+b;
RPP = (a-b)./DL;
TPP = -2*ro1*sqrt(1-(p*VP1).^2)/VP1.*(1-2*(p*VS2).^2)./DL;
TPS = -4*ro1*sqrt(1-(p*VP1).^2)/VP1*VS2^2.*p.*sqrt(1-(p*VP2).^2)/VP2./DL;

figure(1)
subplot(2,2,1), hold on
plot(th1*deg,abs(RPP),'r','Linewidth',3)
plot(th1*deg,abs(TPP),'m-.',th1*deg,abs(TPS),'g:','Linewidth',3)
plot(th1c*deg,0,'ro','Linewidth',3)
plot(th1*deg,abs(RPP),'r','Linewidth',3)
title('COEF. REFLEXION-TRANSMISSION LIQUIDE-SOLIDE'), xlabel('THETA1 (DEG)'), ylabel('MODULE'), 
legend('EAU - SEDIMENT MOU',2)
axis([0,90,0,1.2]), grid on, box on, set(gca,'Color',[1. 0.8 0.5])

subplot(2,2,3), hold on
plot(th1*deg,angle(RPP)*deg,'r','Linewidth',3)
plot(th1*deg,angle(TPP)*deg,'m-.',th1*deg,angle(TPS)*deg,'g:','Linewidth',3)
plot(th1c*deg,-200,'ro','Linewidth',3)
title(['VP2 = ' num2str(VP2) ' M/S , VS2 = ' num2str(VS2) ' M/S , VP2/VS2 = ' num2str(VP2/VS2) ' , RHO2 = ' num2str(ro2)]), xlabel('THETA1 (DEG)'), ylabel('PHASE (DEG)')
legend('RPP','TPP','TPS','thetac(VP1/VP2)',2) 
axis([0,90,-200,200]), grid on, box on, set(gca,'Color',[1. 0.8 0.5])

% cas sediment consolide
VP2 = 2800; VS2 = VP2/1.7; 
ro1 = 1; ro2 = 2;

th1c = asin(VP1/VP2); th2c = asin(VP1/VS2);

p = sin(th1)/VP1;

a = -ro2*sqrt(1-(p*VP1).^2)/VP1.*((1-2*(p*VS2).^2).^2+4*VS2^4*p.^2.*sqrt(1-(p*VP2).^2)/VP2.*sqrt(1-(p*VS2).^2)/VS2);
b = -ro1*sqrt(1-(p*VP2).^2)/VP2;

DL = a+b;
RPP = (a-b)./DL;
TPP = -2*ro1*sqrt(1-(p*VP1).^2)/VP1.*(1-2*(p*VS2).^2)./DL;
TPS = -4*ro1*sqrt(1-(p*VP1).^2)/VP1*VS2^2.*p.*sqrt(1-(p*VP2).^2)/VP2./DL;

figure(1)
subplot(2,2,2), hold on
plot(th1*deg,abs(RPP),'r','Linewidth',3)
plot(th1*deg,abs(TPP),'m-.',th1*deg,abs(TPS),'g:','Linewidth',3)
plot(th1c*deg,0,'ro',th2c*deg,0,'go','Linewidth',3)
plot(th1*deg,abs(RPP),'r','Linewidth',3)
title(['VP1 = ' num2str(VP1) ' M/S , RHO1 = ' num2str(ro1)]), xlabel('THETA1 (DEG)'), ylabel('MODULE')
legend('EAU - SEDIMENT CONSOLIDE',2)
axis([0,90,0,1.2]), grid on, box on, set(gca,'Color',[1. 0.8 0.5])

subplot(2,2,4), hold on
plot(th1*deg,angle(RPP)*deg,'r','Linewidth',3)
plot(th1*deg,angle(TPP)*deg,'m-.',th1*deg,angle(TPS)*deg,'g:','Linewidth',3)
plot(th1c*deg,-200,'ro',th2c*deg,-200,'go','Linewidth',3)
title(['VP2 = ' num2str(VP2) ' M/S , VS2 = ' num2str(VS2,4) ' M/S , VP2/VS2 = ' num2str(VP2/VS2) ' , RHO2 = ' num2str(ro2)]), xlabel('THETA1 (DEG)'), ylabel('PHASE (DEG)')
legend('RPP','TPP','TPS','thetac(VP1/VP2)','thetac(VP1/VS2)',2) 
axis([0,90,-200,200]), grid on, box on, set(gca,'Color',[1. 0.8 0.5])