%STRATIF calcule ondes harmoniques montante + descendante dans stratification
% coef. reflexion par rapport montante/descendante 
% et trace sismique par TF inverse
% voir : http://jmmeost.free.fr/Sismique/osis.html#CS
% fig  = http://jmmeost.free.fr/Sismique/stratif.gif

clear all, clf

dz = 20; % deltaz pour representation onde
HC =  [   0,1000, 500, 500,1000]; NC = length(HC); ZC = cumsum(HC); % epaisseur couches
VP =  [1500,1500,2000,2300,2500]; % vitesses 2500
RHO = [1000,1000,1800,2000,2200]; % densites 2200
IP = RHO.*VP; %impedances
windowSize = 2; IP12 = 2*filter(ones(1,windowSize)/windowSize,1,IP); %I1+I2, I2+I3...
RP = -diff(IP)./IP12(2:NC); RP(1)=0; RP = [RP,0]; TP = 1+RP; %coef. reflexion

subplot(1,4,1); hold on; grid on; box on;
stairs(RHO,ZC,'b'); stairs(VP,ZC,'r'); stairs(IP*1e-4,ZC,'g'); 
legend('RHO','VP','IP'); title('VITESSES , DENSITE, IMPEDANCE');  xlabel('M/S , KG/M3, KG/M2S *10-4'); ylabel('Z (M)'); 
axis([0 3000 0 3000]); set(gca,'YDir','reverse');

dt = 4e-3; nt = 1024; tt = [0:dt:(nt-1)*dt]; % echantillonnage temps
df = 1/nt/dt; fn = 1/2/dt; ff = 0:df:fn; om = 2*pi*ff; %echantillonnage frequence
fc = 30; fe = 15; % frequence centrale et ecart type de la gaussienne
FG = [0 exp(-((ff(2:fix((fc+2*fe)/df))-fc)/fe).^2)]; % spectre gaussien pour ondelette

iom = 0;
for fn = df:df:(length(FG)-1)*df %frequences utilisees par ondelette
om = 2*pi*fn; iom = iom+1;
KZC = om./VP; KHC = KZC.*HC; % kzj , kzj*hj
EKH1 = exp(i*KHC(NC-1)); EKH2 = exp(i*KHC(NC-2)); EKH3 = exp(i*KHC(NC-3)); % exp(i*kzj*hj)
Z0= 0:dz:HC(NC); Z1= 0:dz:HC(NC-1); Z2= 0:dz:HC(NC-2); Z3= 0:dz:HC(NC-3); % zj
EKZ0 = exp(i*KZC(NC).*Z0); EKZ1 = exp(i*KZC(NC-1).*Z1); EKZ2 = exp(i*KZC(NC-2).*Z2); EKZ3 = exp(i*KZC(NC-3).*Z3); % exp(i*kzj*zj)
% couche 1 et 1/2 espace
D1 = 1;         DZ1 = D1*conj(EKZ1); % descendante
M1 = RP(NC-1);  MZ1 = M1*EKZ1; % montante
R1 = M1/D1; RZ1 = MZ1./DZ1; % coef. reflexion a interface 1/0
D0 = TP(NC-1); DZ0 = D0*EKZ0;  %descendante dans 1/2 espace
% couche 2
D2 = (conj(EKH1)+RP(NC-2)*RP(NC-1)*EKH1)*D1/TP(NC-2); DZ2 = D2*conj(EKZ2);  
M2 = (RP(NC-2)*conj(EKH1)+RP(NC-1)*EKH1)*D1/TP(NC-2); MZ2 = M2*EKZ2; 
RH1 = M2/D2; RHZ1 = MZ2./DZ2.*conj(EKZ2).*conj(EKZ2); %coef. reflexion sur couche 1
R2 =  (RH1*conj(EKH1)-R1*EKH1)/(conj(EKH1)-R1*RH1*EKH1); %coef. reflexion sur interface 2/1
RZ2 = (RHZ1*conj(EKH1)-R1*EKH1)./(conj(EKH1)-R1*RHZ1*EKH1); 
% couche 3
D3 = (conj(EKH2)+RP(NC-3)*RH1*EKH2)*D2/TP(NC-3); DZ3 = D3*conj(EKZ3);
M3 = (RP(NC-3)*conj(EKH2)+RH1*EKH2)*D2/TP(NC-3); MZ3 = M3*EKZ3; 
RH12 = M3/D3; RHZ2 = MZ3./DZ3.*conj(EKZ3).*conj(EKZ3); %coef. reflexion sur couches 2 et 1
R3 = (RH12*conj(EKH2)-RH1*EKH2)/(conj(EKH2)-RH1*RH12*EKH2); %coef. reflexion sur interface 3/2
RZ3 = (RHZ2*conj(EKH2)-RH1*EKH2)./(conj(EKH2)-RH1*RHZ2*EKH2);

if iom == 10
subplot(1,4,2); hold on; grid on; box on;
plot(real(DZ3),ZC(NC-3)-Z3,'r',real(MZ3),ZC(NC-3)-Z3,'b',real(DZ3+MZ3),ZC(NC-3)-Z3,'g');
plot(real(DZ2),ZC(NC-2)-Z2,'r',real(MZ2),ZC(NC-2)-Z2,'b',real(DZ2+MZ2),ZC(NC-2)-Z2,'g');
plot(real(DZ1),ZC(NC-1)-Z1,'r',real(MZ1),ZC(NC-1)-Z1,'b',real(DZ1+MZ1),ZC(NC-1)-Z1,'g');
plot(real(DZ0),ZC(NC-1)+Z0,'g');
legend(['D (f =' num2str(fn) 'HZ)'],'M','D+M'); title('DESC. +  MONTANTE');
axis([-10 10 0 3000]); set(gca,'YDir','reverse');

subplot(1,4,3); hold on; grid on; box on;
plot(RP,ZC,'or');
plot(-abs(RHZ2),ZC(NC-3)-Z3,'g', -abs(RZ3),ZC(NC-3)-Z3,'b', real(R3),ZC(NC-3),'*b');
plot(-abs(RHZ1),ZC(NC-2)-Z2,'g', -abs(RZ2),ZC(NC-2)-Z2,'b', real(R2),ZC(NC-2),'*b');
plot( -abs(RZ1),ZC(NC-1)-Z1,'g',  real(R1),ZC(NC-1),'*b');
legend('R VRAI','|M/D|','R(Z)','R'); title('COEF. REFLEXION PAR M/D'); xlabel('');
axis([-1 0 0 3000]); set(gca,'YDir','reverse');
end

FF(iom+1) = RH12*EKH3*EKH3*FG(iom+1); % onde reflechie par couches 2 et 1 prolongee en z = 0
FF1(iom+1) = RH12*EKH3*EKH3*FF(iom+1); % réflexion de la precedente sur surface libre (R = 1)
FF2(iom+1) = RH12*EKH3*EKH3*FF1(iom+1); % réflexion de la precedente sur surface libre (R = 1)
end

ZE = zeros(nt/2+1-length(FF)); FP = [FF,ZE(1,:)]; FP1 = [FF1,ZE(1,:)]; FP2 = [FF2,ZE(1,:)];% frequences non utilisees pour ondelette
FN = conj(fliplr(FP));  FN1 = conj(fliplr(FP1)); FN2 = conj(fliplr(FP2));% frequences negatives pour trace reelle
TF = [FP,FN(2:nt/2)]; TF1 = [FP1,FN1(2:nt/2)]; TF2 = [FP2,FN2(2:nt/2)];% toutes les frequences

TR = ifft(conj(TF)); TR1 = ifft(conj(TF1)); TR2 = ifft(conj(TF2));% TF inverse
TR = TR/max(abs(TR)); TR1 = TR1/max(abs(TR)); TR2 = TR2/max(abs(TR)); % normalisation par max(TR)

subplot(1,4,4); hold on; grid on; box on;
it1 = fix(4*HC(2)/VP(2)/dt); nt1 = [it1-10:nt]; it2 = fix(6*HC(2)/VP(2)/dt); nt2 = [it2-10:nt]; % indices pour reflexions surface
plot(5*real(TR),tt,'b',5*real(TR1(nt1)),tt(nt1),'r',5*real(TR2(nt2)),tt(nt2),'g');
ylabel('T (S)'); title('TRACE SISMIQUE'); legend('R. PAR STRATIF.','R. 1 FOIS PAR Z=0','R. 2 FOIS PAR Z=0',['|D| = exp((F-' int2str(fc) '/' int2str(fe) ')**2)']);
axis([-1 1 0 nt*dt]); set(gca,'YDir','reverse');