%TRANSTRAT calcule les premiers termes du coef. de transmission d'une stratification cyclique
% voir : http://jmmeost.free.fr/Sismique/osis.html#SC
% fig1 = http://jmmeost.free.fr/Sismique/transtrat.gif  

clear all, clf
  
RE = 0.1; TE=1+RE; RS = -RE; TS = 1-RE; % coef. sur interface

M = 500; dm = 10;
for m = 2:dm:M %2m+1 couches 2m+2 interfaces

N1 = zeros(2*m+1); N1(1) = 1 ; N2 = zeros(2*m+1); N2(1) = 1 ; 
N3 = zeros(2*m+1); N3(1) = 1 ; N4 = zeros(2*m+1); N4(1) = 1 ;

for n = 2:2*m+1
   N1(n) = n; % trajets avec 2 reflexions
   N2(n) = N1(n)+N2(n-1); % trajets avec 4 reflexions
   N3(n) = N2(n)+N3(n-1); % trajets avec 6 reflexions
   N4(n) = N3(n)+N4(n-1); % trajets avec 8 reflexions
end

T(m) = (1-RE^2)^(m+1); %transmission directe
T1(m) = T(m).*N1(2*m+1)*RE^2; % 1 trajet multiple 
T2(m) = T(m).*(N2(2*m+1)*RE^4-N1(2*m)*RE^2*(1-RE^2)); % 2 trajets multiples
T3(m) = T(m).*(N3(2*m+1)*RE^6-(2*N2(2*m+1)-2)*(1-RE^2)*RE^4+N1(2*m-1)*((1-RE^2)^2)*RE^2); % 3 trajets multiples
T4(m) = T(m).*(N4(2*m+1)*RE^8-(3*N3(2*m+1)+N2(2*m)-3)*(1-RE^2)*RE^6+(3*N2(2*m+1)-2*N2(2)-2)*(1-RE^2)^2*RE^4-N1(2*m-2)*(1-RE^2)^3*RE^2); % 4 trajets multiples

end

m = 2:dm:M;
plot(2*m+1,T(m),'b',2*m+1,T1(m),'r',2*m+1,T2(m),'g',2*m+1,T3(m),'m',2*m+1,T4(m),'c')
legend('TS0 R=0.1','TS1 ','TS2','TS3','TS4')
title('COEFFICIENT TRANSMISSION ONDE PLANE SH A TRAVERS STRATIFICATION CYCLIQUE') 
xlabel('NOMBRE COUCHES'), ylabel('MODULE TS')