%REFLSH calcule coef. reflexion onde SH en fonction angle incidence
% et montre la transition transmission/evanescence a l'incidence critique
% voir : http://jmmeost.free.fr/Sismique/osis.html#SH
% fig1 = http://jmmeost.free.fr/Sismique/reflsh.gif.
% fig2 = http://jmmeost.free.fr/Sismique/reflsht.gif

clear all, clf

V1 = 2000 ; V2 = 2600; %vitesses
eta1 = 0:.001:pi/2; deg = 180/pi; %incidences
p = sin(eta1)/V1; etac = asin(V1/V2)*deg %parametre, incidence critique

% coefficients
I1 = V1*sqrt(1-(p*V1).^2); I2 = V2*sqrt(1-(p*V2).^2); %impedance (rho constant)
R = (I1-I2)./(I1+I2); T = 2*I1./(I1+I2); %coef. reflexion-transmission

figure(1)
subplot(2,1,1), 
plot(eta1*deg,abs(T),'r',eta1*deg,abs(R),'b')
title('REFLEXION ONDE SH V1 = 2000 M/S V2 = 2600 M/S RHO = cte')
ylabel('MODULE'), legend('T','R')
subplot(2,1,2), 
plot(eta1*deg,angle(T)*deg,'r',eta1*deg,angle(R)*deg,'b')
xlabel('ETA1 (DEG)'), ylabel('PHASE (DEG)')

% ondes
om = 2*pi*2; ip = 0 % f = 2 Hz
for eta1 = [45,50,80,52] %incidences pre et post critique
ip = ip+1;
p = sin(eta1*pi/180)/V1; 
I1 = V1*sqrt(1-(p*V1).^2); I2 = V2*sqrt(1-(p*V2).^2);
R = (I1-I2)./(I1+I2); T = 2*I1./(I1+I2);

dx = 50; dz = 50; %delta x et z en m
x = 0:dx:6000; z = -3000:dz:0; [X1,Z1] = meshgrid(x,z);
px = p*X1; pz1 = Z1*I1/V1^2;
uyi = exp(i*om*(px+pz1)); uyr = R*exp(i*om*(px-pz1)); %ondes incidente et reflechie

z = 0:dz:3000; [X2,Z2] = meshgrid(x,z);
pz2 = Z2*I2/V2^2;
uyt = T*exp(i*om*(px+pz2)); %onde transmise

figure(2)

subplot(2,2,ip), hold on
pcolor([X1;X2],[Z1;Z2],real([uyi+uyr;uyt])), plot([0,6000],[0,0],'k')
set(gca, 'YDir', 'reverse','DataAspectRatio',[1,1,1]),shading('interp')
switch ip
case 1
   title('REFLEXION PARTIELLE ETA1 = 45°'),ylabel('Z (M)')
case 2
   title('ETA1 = 50° < ETAC = 50.3°')
   ylabel('TRANSMISE | INCID. + REFL.')
case 3
   title('REFLEXION TOTALE ETA1 = 80°'), xlabel('X (M)'),ylabel('Z (M)')
case 4
   title('REFLEXION TOTALE ETA1 = 52°'), xlabel('X (M)')
   ylabel('EVANESCENTE | STATIONNAIRE')
   gtext('REFLEXION-TRANSMISSION ONDE PLANE HARMONIQUE SH EN Z = 0 V1 = 2000 M/S V2 = 2600 M/S  F = 2 HZ')
end
end