%MONITORPSV calcule ondes harmoniques montante + descendante P et SV dans stratification  
% et trace sismique par TF inverse pour deux modeles A et B
% voir : http://jmmeost.free.fr/Sismique/osis.html#MS
% fig1  = http://jmmeost.free.fr/Sismique/monitormod.gif
% fig2  = http://jmmeost.free.fr/Sismique/monitorrefl.gif
% fig3  = http://jmmeost.free.fr/Sismique/monitorpsv0.gif
% fig4  = http://jmmeost.free.fr/Sismique/monitorpsv1.gif
% fig5  = http://jmmeost.free.fr/Sismique/monitorpsv2.gif

clear all

% NC-1 couches de 2 (eau) a NC (1/2 espace inf) ; NC-1 interfaces de 1 (surface libre) a NC-1 
HC =  [    0,    300,    200,     500,      40,      350,   1000]; NC = length(HC); ZC = cumsum(HC); 

% modele vitesses et densite  A et B
VPA =  [ 1500,   1500,   2500,    2900,     3200,    3000,   3600];
VPB =  [ 1500,   1500,   2500,    2900,     2800,    3000,   3600];

VSA =  [    0,      0,   1210,   1700,      1550,    1800,   2200];
VSB =  [    0,      0,   1210,   1700,      1550,    1800,   2200];

RHOA = [ 1000,   1000,   2400,   2500,      2300,    2200,   2400];
RHOB = [ 1000,   1000,   2400,   2500,      2300,    2200,   2400];

% spectre pour reponse impulsionnelle
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

VP = VPA; VS = VSA; RHO = RHOA; VPVSA = VP./VS;  IPA = RHO.*VP; ISA = RHO.*VS;
IP = IPA; IP12 = 2*filter(ones(1,2)/2,1,IP); % impedances P, IP12 = IP1+IP2...
RP0A = -diff(IP)./IP12(2:NC); RP0A(1)=1; TP0A = 2*IP(2:NC-1)./IP12(3:NC); % coef. refl. trans. P incidence verticale
IS = ISA; IS12 = 2*filter(ones(1,2)/2,1,IS); % impedances S, IS12 = IS1+IS2...
RS0A = -diff(IS)./IS12(2:NC); RS0A(2)=1; TS0A = 2*IS(2:NC-1)./IS12(3:NC); % coef. refl. trans. S incidence verticale

VP = VPB; VS = VSB; RHO = RHOB; VPVSB = VP./VS; IPB = RHO.*VP; ISB = RHO.*VS;
IP = IPB; IP12 = 2*filter(ones(1,2)/2,1,IP); % impedances P, IP12 = IP1+IP2...
RP0B = -diff(IP)./IP12(2:NC); RP0B(1)=1; TP0B = 2*IP(2:NC-1)./IP12(3:NC); % coef. refl. trans. P incidence verticale
IS = ISB; IS12 = 2*filter(ones(1,2)/2,1,IS); % impedances S, IS12 = IS1+IS2...
RS0B = -diff(IS)./IS12(2:NC); RS0B(2)=1; TS0B = 2*IS(2:NC-1)./IS12(3:NC); % coef. refl. trans. S incidence verticale

figure(1), clf
subplot(1,2,1); hold on; grid on; box on;
stairs(RHOA,ZC,'g'); stairs(VPA,ZC,'r'); stairs(VSA,ZC,'b');  stairs(1000*(VPVSA-1.7),ZC,'m');
stairs(RHOB,ZC,'g:'); stairs(VPB,ZC,'r:'); stairs(VSB,ZC,'b:');  stairs(1000*(VPVSB-1.7),ZC,'m:');
legend('RHOA','VPA','VSA','(VPA/VSA-1.7)*1000','RHOB','VPB','VSB','(VPB/VSB-1.7)*1000'); 
xlabel('VITESSES (M/S) , DENSITE (KG/M3)');  ylabel('PROFONDEUR (M)'); title('MONITORING SISMIQUE : MODELES A ET B'); 
axis([0 3500 0 1500]); set(gca,'YDir','reverse');

subplot(1,2,2); hold on; grid on; box on;
stairs(IPA*1e-7,ZC,'r'); stairs(ISA*1e-7,ZC,'b');
plot(RP0A,ZC(1:NC-1),'^r','MarkerFaceColor','r'); plot(RS0A(3:NC-1),ZC(3:NC-1),'^b','MarkerFaceColor','b');
plot(TP0A,ZC(2:NC-1),'vr','MarkerFaceColor','r'); plot(TS0A,ZC(2:NC-1),'vb','MarkerFaceColor','b');
stairs(IPB*1e-7,ZC,'r:'); stairs(ISB*1e-7,ZC,'b:');
plot(RP0B,ZC(1:NC-1),'^r','MarkerSize',12); plot(RS0B(3:NC-1),ZC(3:NC-1),'^b','MarkerSize',12);
plot(TP0B,ZC(2:NC-1),'vr','MarkerSize',12); plot(TS0B,ZC(2:NC-1),'vb','MarkerSize',12);
legend('IPA*1e-7','ISA*1e-7','RPPA','RSSA','TPPA','TSSA','IPB*1e-7','ISB*1e-7','RPPB','RSSB','TPPB','TSSB'); 
xlabel('COEF. REFLEXION , IMPEDANCES (KG/M2S) , COEF. TRANSMISSION'); title('INCIDENCE VERTICALE');
axis([-.4 1.2 0 1500]); set(gca,'YDir','reverse');

% modele A

VP = VPA; VS = VSA; RHO = RHOA;

figure(2), clf, figure(3), clf, figure(4), clf, figure(5), clf

% angle incidence dans l'eau
ith = 0;
for thd = [0 15 22]
ith = ith+1;
th = thd*pi/180; p = sin(th)/VP(2);

% calcul des coef. reflexion-transmission a interface ic / ic+1

% interface eau (ic = 2) - sediment (ic = 3)
ic = 2;
VP1 = VP(ic); VP2 = VP(ic+1); VS2 = VS(ic+1); ro1 = RHO(ic); ro2 = RHO(ic+1);
qP1 = sqrt(1-(p*VP1).^2)/VP1; qP2 = sqrt(1-(p*VP2).^2)/VP2; qS2 = sqrt(1-(p*VS2).^2)/VS2; % cos(theta)/V
cS2 = 1-2*(p*VS2).^2; % cos(2theta)
qPH(ic) = qP1*HC(ic); qPH(ic+1) = qP2*HC(ic+1); qSH(ic+1) = qS2*HC(ic+1); % hcos(theta)/V

A1 = [-qP1 ; ro1 ; 0]; A2 = [-qP2 ; -ro2*cS2 ; 2*ro2*VS2^2*p.*qP2]; A3 = [-p ; -2*ro2*VS2^2*p.*qS2 ; -ro2*cS2];
A = [A1, A2, A3]; 

RP = A\[-qP1; -ro1 ; 0]; % onde incidente P descendante
RPP(ic) = RP(1); RPS(ic) = 0; TPP(ic) = RP(2); TPS(ic) = RP(3);

rP = A\[-qP2; ro2*cS2; 2*ro2*VS2^2*p.*qP2]; % onde incidente P montante
tPP(ic) = rP(1); tPS(ic) = 0; rPP(ic) = rP(2); rPS(ic) = rP(3);

RSP(ic) = 0; RSS(ic) = 0; TSP(ic) = 0; TSS(ic) = 0; % pas d'onde S descendante

rS = A\[p ; -2*ro2*VS2^2*p.*qS2 ; ro2*cS2]; % onde incidente S montante
tSP(ic) = rS(1);  tSS(ic) = 0; rSP(ic) = rS(2); rSS(ic) = rS(3);

% stratification elastique

for ic = 3:NC-1

VP1 = VP(ic); VP2 = VP(ic+1); VS1 = VS(ic); VS2 = VS(ic+1); ro1 = RHO(ic); ro2 = RHO(ic+1);
qP1 = sqrt(1-(p*VP1).^2)/VP1; qS1 = sqrt(1-(p*VS1).^2)/VS1; qP2 = sqrt(1-(p*VP2).^2)/VP2; qS2 = sqrt(1-(p*VS2).^2)/VS2; %cos(theta)/V
cS1 = 1-2*(p*VS1).^2; cS2 = 1-2*(p*VS2).^2; %cos(2theta)
qPH(ic+1) = qP2*HC(ic+1); qSH(ic+1) = qS2*HC(ic+1); % hcos(theta)/V
    
A1 = [p ; -qP1 ; ro1*cS1 ; 2*ro1*VS1^2*p.*qP1];   A2 = [qS1 ; p ; -2*ro1*VS1^2*p.*qS1 ; ro1*cS1];
A3 = [-p ; -qP2 ; -ro2*cS2 ; 2*ro2*VS2^2*p.*qP2]; A4 = [qS2 ; -p ; -2*ro2*VS2^2*p.*qS2 ; -ro2*cS2];
A = [A1 , A2, A3, A4]; 

RP = A\[-p; -qP1; -ro1*cS1 ; 2*ro1*VS1^2*p.*qP1]; % onde incidente P descendante
RPP(ic) = RP(1); RPS(ic) = RP(2); TPP(ic) = RP(3); TPS(ic) = RP(4);

rP = A\[p; -qP2; ro2*cS2 ; 2*ro2*VS2^2*p.*qP2]; % onde incidente P montante
tPP(ic) = rP(1); tPS(ic) = rP(2); rPP(ic) = rP(3); rPS(ic) = rP(4);

RS = A\[qS1 ; -p ; -2*ro1*VS1^2*p.*qS1 ; -ro1*cS1]; % onde incidente S descendante
RSP(ic) = RS(1); RSS(ic) = RS(2); TSP(ic) = RS(3); TSS(ic) = RS(4);

rS = A\[qS2 ; p ; -2*ro2*VS2^2*p.*qS2 ; ro2*cS2]; % onde incidente S montante
tSP(ic) = rS(1); tSS(ic) = rS(2); rSP(ic) = rS(3); rSS(ic) = rS(4);

end

TTPP = 2*cumsum(qPH(3:NC-1)); % temps trajet vertical P 2Sum(Hcostheta/VP)
TTSS = 2*cumsum(qSH(3:NC-1)); % temps trajet vertical S 2Sum(Hcoseta/VS)
TTPS = (TTPP+TTSS)/2;

figure(2)
switch thd
case 0
subplot(2,3,1), hold on; grid on; box on; ylabel('PROFONDEUR (M)'); title(['INCIDENCE ' int2str(thd) '°']);
case 15
subplot(2,3,2), hold on; grid on; box on; title(['INCIDENCE ' int2str(thd) '°']);
case 22
subplot(2,3,3), hold on; grid on; box on; title(['INCIDENCE ' int2str(thd) '°']);
end

plot(RPP(2:NC-1),ZC(2:NC-1),'^r','MarkerFaceColor','r'), plot(RSS(2:NC-1),ZC(2:NC-1),'^b','MarkerFaceColor','b')
plot(RPS(2:NC-1),ZC(2:NC-1),'^g','MarkerFaceColor','g'), plot(RSP(2:NC-1),ZC(2:NC-1),'^m','MarkerFaceColor','m')
if thd == 0
plot(RPP(2),ZC(2),'^r','MarkerSize',12), plot(RSS(2),ZC(2),'^b','MarkerSize',12), 
plot(RPS(2),ZC(2),'^g','MarkerSize',12), plot(RSP(2),ZC(2),'^m','MarkerSize',12),
legend('RPPA','RSSA','RPSA','RSPA','RPPB','RSSB','RPSB','RSPB')
end
axis([-0.3 0.3 0 1500]); set(gca,'YDir','reverse'); xlabel('COEF. REFLEXION');

switch thd
case 0
subplot(2,3,4), hold on; grid on; box on; ylabel('PROFONDEUR (M)'); 
case 15
subplot(2,3,5), hold on; grid on; box on; 
case 22
subplot(2,3,6), hold on; grid on; box on; 
end

plot(TPP(2:NC-1),ZC(2:NC-1),'vr','MarkerFaceColor','r'), plot(TSS(2:NC-1),ZC(2:NC-1),'vb','MarkerFaceColor','b')
plot(TPS(2:NC-1),ZC(2:NC-1),'vg','MarkerFaceColor','g'), plot(TSP(2:NC-1),ZC(2:NC-1),'vm','MarkerFaceColor','m')
if thd == 0
plot(TPP(2),ZC(2),'vr','MarkerSize',12), plot(TSS(2),ZC(2),'vb','MarkerSize',12), 
plot(TPS(2),ZC(2),'vg','MarkerSize',12), plot(TSP(2),ZC(2),'vm','MarkerSize',12),
legend('TPPA','TSSA','TPSA','TSPA','TPPB','TSSB','TPSB','TSPB')
end
axis([-0.5 1.2 0 1500]); set(gca,'YDir','reverse'); xlabel('COEF. TRANSMISSION');

% calcul pour frequences utilisees par ondelette
iom = 0;
for fo = df:df:(length(FG)-1)*df 
om = 2*pi*fo; iom = iom+1;

DM = eye(4); % matrice de transfert descendantes P et S - montantes P et S

EKPH = exp(i*om*qPH(2:NC)); CKPH = conj(EKPH); % propagation P dans couches
EKSH = exp(i*om*qSH(2:NC)); CKSH = conj(EKSH); % propagation S dans couches

% matrice transfert a travers couches et stratification

for ic = 3:NC-2
    
div = TPP(ic)*TSS(ic)-TPS(ic)*TSP(ic);
DP = [TSS(ic)*CKPH(ic) , -TSP(ic)*CKSH(ic) , (TSP(ic)*rPS(ic)-TSS(ic)*rPP(ic))*EKPH(ic) , (TSP(ic)*rSS(ic)-TSS(ic)*rSP(ic))*EKSH(ic)]/div;
DS = [TPS(ic)*CKPH(ic) , -TPP(ic)*CKSH(ic) , (TPP(ic)*rPS(ic)-TPS(ic)*rPP(ic))*EKPH(ic) , (TPP(ic)*rSS(ic)-TPS(ic)*rSP(ic))*EKSH(ic)]/(-div);
MP = RPP(ic)*DP + RSP(ic)*DS + [0 , 0 , tPP(ic)*EKPH(ic) , tSP(ic)*EKSH(ic)];
MS = RPS(ic)*DP + RSS(ic)*DS + [0 , 0 , tPS(ic)*EKPH(ic) , tSS(ic)*EKSH(ic)];

DM = DM*[DP ; DS ; MP ; MS];
end

DM = DM*[1 , 0 ; 0 , 1 ; RPP(NC-1) , RSP(NC-1) ; RPS(NC-1) , RSS(NC-1)]; % interface avec 1/2 espace inferieur

% cas sismique reflexion terrestre
% onde incidente P d'amplitude 1 au sommet stratification
DP1 = DM(1:2,:)\[1 ; 0]; %DP1 dans couche du bas 
RHPPS = DM(3:4,:)*DP1;  % ondes montantes P et S au sommet de la stratification
RHPP = RHPPS(1); RHPS = RHPPS(2); %coef. reflexion sur stratification elastique

AA(iom+1) =  RHPP*EKPH(2)*EKPH(2)*FG(iom+1); % ondes P sous fond de l'eau
BB(iom+1) =  RHPS*EKPH(2)*EKSH(2)*FG(iom+1); % ondes S sous fond de l'eau

% onde incidente S d'amplitude 1 au sommet stratification
DS1 = DM(1:2,:)\[0 ; 1]; %DP1 dans couche du bas 
RHSPS = DM(3:4,:)*DS1;  % ondes montantes P et S au sommet de la stratification
RHSP = RHSPS(1); RHSS = RHSPS(2); %coef. reflexion sur stratification elastique

CC(iom+1) =  RHSP*EKSH(2)*EKPH(2)*FG(iom+1); % ondes P sous fond de l'eau
DD(iom+1) =  RHSS*EKSH(2)*EKSH(2)*FG(iom+1); % ondes S sous fond de l'eau

% cas sismique reflexion marine
% matrice transfert a travers couche et stratification sous eau
DP = [CKPH(2) , 0 , -rPP(2)*EKPH(2) , -rSP(2)*EKSH(2)]/TPP(2);
MP = RPP(2)*DP + [0 , 0, tPP(2)*EKPH(2) , tSP(2)*EKSH(2)];
DME = [DP ; MP]*DM;

DE1 = DME(1,:)\1; %DP1 et DS1 pour onde incidente P d'amplitude 1 dans l'eau
RHEE = DME(2,:)*DE1; %coef. reflexion sur stratification elastique sous eau
RHEPS = DM(3:4,:)*DE1; RHEP = RHEPS(1); RHES = RHEPS(2); % ondes montantes P et S dans couche sous eau

FF(iom+1) =  RHEE*FG(iom+1); % onde P au fond de l'eau reflechie par stratification 
GG(iom+1) =  RHEP*EKPH(2)*FG(iom+1); % ondes P sous fond de l'eau reflechie par stratification
HH(iom+1) =  RHES*EKSH(2)*FG(iom+1); % ondes S sous fond de l'eau reflechie par stratification

end % frequences

ZE = zeros(nt/2+1-length(FF)); 

FP = [AA,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTPPA(ith) = max(abs(TF)); TPPA(ith,:) = TF/MTPPA(ith);
FP = [BB,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTPSA(ith) = max(abs(TF)); TPSA(ith,:) = TF/MTPSA(ith);
FP = [CC,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTSPA(ith) = max(abs(TF)); TSPA(ith,:) = TF/MTSPA(ith);
FP = [DD,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTSSA(ith) = max(abs(TF)); TSSA(ith,:) = TF/MTSSA(ith);
FP = [FF,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTEEA(ith) = max(abs(TF)); TEEA(ith,:) = TF/MTEEA(ith);
FP = [GG,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTEPA(ith) = max(abs(TF)); TEPA(ith,:) = TF/MTEPA(ith);
FP = [HH,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTESA(ith) = max(abs(TF)); TESA(ith,:) = TF/MTESA(ith);

switch thd
    case 0, figure(3)
    case 15, figure(4)
    case 22, figure (5)
end

subplot(1,5,1), hold on; grid on; box on; 
for ic = 1:NC-3  plot([-1,-0.5],[TTPP(ic),TTPP(ic)],'r'); plot([0.5,1],[TTPS(ic)-TTPS(1)+TTPP(1),TTPS(ic)-TTPS(1)+TTPP(1)],'g'); end
for ic = 2:NC-3  plot([1,1.5],[TTPS(ic)-TTPS(2)+TTPP(2),TTPS(ic)-TTPS(2)+TTPP(2)],'g'); end 
plot(TPPA(ith,:),tt,'r'); 
xlabel('P-P SANS EAU'); legend('TTPP','TTPP+PS');
ylabel('TEMPS DE PROPAGATION VERTICALE DE L"ONDE PLANE DANS STRATIFICATION (S)');

subplot(1,5,2), hold on; grid on; box on; 
for ic = 1:NC-3  plot([-1,-0.5],[TTSS(ic),TTSS(ic)],'b'); plot([0.5,1],[TTPS(ic)-TTPS(1)+TTSS(1),TTPS(ic)-TTPS(1)+TTSS(1)],'g'); end
for ic = 2:NC-3 plot([1,1.5],[TTPS(ic)-TTPS(2)+TTSS(2),TTPS(ic)-TTPS(2)+TTSS(2)],'g'); end
plot(TSSA(ith,:),tt,'b'); 
xlabel('SV-SV SANS EAU'); legend('TTSS','TTSS+PS');

subplot(1,5,3), hold on; grid on; box on;
plot(10*TEEA(ith,:),tt,'k');
xlabel('P-P DANS EAU');

subplot(1,5,4), hold on; grid on; box on;
for ic = 1:NC-3  
plot([-1,-0.5],[TTPP(ic),TTPP(ic)],'r'); plot([0.5,1],[TTPP(1)+TTPP(ic),TTPP(1)+TTPP(ic)],'r'); 
plot([1,1.5],[TTPP(2)+TTPP(ic),TTPP(2)+TTPP(ic)],'r'); plot([2.5,3],[TTPP(4)+TTPP(ic),TTPP(4)+TTPP(ic)],'r');
end
plot(TEPA(ith,:),tt,'r'); 
xlabel('P DANS - P SOUS EAU'); legend('TTPP','TTPP+PP'); %legend(['MAX = ' num2str(100*MTEP,'%3.2g')]);

subplot(1,5,5), hold on; grid on; box on;
for ic = 1:NC-3  plot([-1,-0.5],[TTPS(ic),TTPS(ic)],'g'); plot([0.5,1],[TTPP(ic)-TTPP(1)+TTPS(1),TTPP(ic)-TTPP(1)+TTPS(1)],'r'); end
for ic = 2:NC-3  plot([1,1.5],[TTPP(ic)-TTPP(2)+TTPS(2),TTPP(ic)-TTPP(2)+TTPS(2)],'r'); end
plot(TESA(ith,:),tt,'g');
xlabel('P DANS - SV SOUS EAU'); legend('TTPS','TTPS+PP'); %legend(['MAX = ' num2str(100*MTES,'%3.2g')]);

end % incidence

% modele B

VP = VPB; VS = VSB; RHO = RHOB;

figure(2), figure(3), figure(4), figure(5),

% angle incidence dans l'eau
ith = 0;
for thd = [0 15 22]
ith = ith+1;
th = thd*pi/180; p = sin(th)/VP(2);

% calcul des coef. reflexion-transmission a interface ic / ic+1

% interface eau (ic = 2) - sediment (ic = 3)
ic = 2;
VP1 = VP(ic); VP2 = VP(ic+1); VS2 = VS(ic+1); ro1 = RHO(ic); ro2 = RHO(ic+1);
qP1 = sqrt(1-(p*VP1).^2)/VP1; qP2 = sqrt(1-(p*VP2).^2)/VP2; qS2 = sqrt(1-(p*VS2).^2)/VS2; % cos(theta)/V
cS2 = 1-2*(p*VS2).^2; % cos(2theta)
qPH(ic) = qP1*HC(ic); qPH(ic+1) = qP2*HC(ic+1); qSH(ic+1) = qS2*HC(ic+1); % hcos(theta)/V

A1 = [-qP1 ; ro1 ; 0]; A2 = [-qP2 ; -ro2*cS2 ; 2*ro2*VS2^2*p.*qP2]; A3 = [-p ; -2*ro2*VS2^2*p.*qS2 ; -ro2*cS2];
A = [A1, A2, A3]; 

RP = A\[-qP1; -ro1 ; 0]; % onde incidente P descendante
RPP(ic) = RP(1); RPS(ic) = 0; TPP(ic) = RP(2); TPS(ic) = RP(3);

rP = A\[-qP2; ro2*cS2; 2*ro2*VS2^2*p.*qP2]; % onde incidente P montante
tPP(ic) = rP(1); tPS(ic) = 0; rPP(ic) = rP(2); rPS(ic) = rP(3);

RSP(ic) = 0; RSS(ic) = 0; TSP(ic) = 0; TSS(ic) = 0; % pas d'onde S descendante

rS = A\[p ; -2*ro2*VS2^2*p.*qS2 ; ro2*cS2]; % onde incidente S montante
tSP(ic) = rS(1);  tSS(ic) = 0; rSP(ic) = rS(2); rSS(ic) = rS(3);

% stratification elastique

for ic = 3:NC-1

VP1 = VP(ic); VP2 = VP(ic+1); VS1 = VS(ic); VS2 = VS(ic+1); ro1 = RHO(ic); ro2 = RHO(ic+1);
qP1 = sqrt(1-(p*VP1).^2)/VP1; qS1 = sqrt(1-(p*VS1).^2)/VS1; qP2 = sqrt(1-(p*VP2).^2)/VP2; qS2 = sqrt(1-(p*VS2).^2)/VS2; %cos(theta)/V
cS1 = 1-2*(p*VS1).^2; cS2 = 1-2*(p*VS2).^2; %cos(2theta)
qPH(ic+1) = qP2*HC(ic+1); qSH(ic+1) = qS2*HC(ic+1); % hcos(theta)/V
    
A1 = [p ; -qP1 ; ro1*cS1 ; 2*ro1*VS1^2*p.*qP1];   A2 = [qS1 ; p ; -2*ro1*VS1^2*p.*qS1 ; ro1*cS1];
A3 = [-p ; -qP2 ; -ro2*cS2 ; 2*ro2*VS2^2*p.*qP2]; A4 = [qS2 ; -p ; -2*ro2*VS2^2*p.*qS2 ; -ro2*cS2];
A = [A1 , A2, A3, A4]; 

RP = A\[-p; -qP1; -ro1*cS1 ; 2*ro1*VS1^2*p.*qP1]; % onde incidente P descendante
RPP(ic) = RP(1); RPS(ic) = RP(2); TPP(ic) = RP(3); TPS(ic) = RP(4);

rP = A\[p; -qP2; ro2*cS2 ; 2*ro2*VS2^2*p.*qP2]; % onde incidente P montante
tPP(ic) = rP(1); tPS(ic) = rP(2); rPP(ic) = rP(3); rPS(ic) = rP(4);

RS = A\[qS1 ; -p ; -2*ro1*VS1^2*p.*qS1 ; -ro1*cS1]; % onde incidente S descendante
RSP(ic) = RS(1); RSS(ic) = RS(2); TSP(ic) = RS(3); TSS(ic) = RS(4);

rS = A\[qS2 ; p ; -2*ro2*VS2^2*p.*qS2 ; ro2*cS2]; % onde incidente S montante
tSP(ic) = rS(1); tSS(ic) = rS(2); rSP(ic) = rS(3); rSS(ic) = rS(4);

end

TTPP = 2*cumsum(qPH(3:NC-1)); % temps trajet vertical P 2Sum(Hcostheta/VP)
TTSS = 2*cumsum(qSH(3:NC-1)); % temps trajet vertical S 2Sum(Hcoseta/VS)
TTPS = (TTPP+TTSS)/2;

figure(2)
switch thd
case 0
subplot(2,3,1), hold on; grid on; box on; ylabel('PROFONDEUR (M)'); title(['INCIDENCE ' int2str(thd) '°']);
case 15
subplot(2,3,2), hold on; grid on; box on; title(['INCIDENCE ' int2str(thd) '°']);
case 22
subplot(2,3,3), hold on; grid on; box on; title(['INCIDENCE ' int2str(thd) '°']);
end

plot(RPP(2:NC-1),ZC(2:NC-1),'^r','MarkerSize',12), plot(RSS(2:NC-1),ZC(2:NC-1),'^b','MarkerSize',12)
plot(RPS(2:NC-1),ZC(2:NC-1),'^g','MarkerSize',12), plot(RSP(2:NC-1),ZC(2:NC-1),'^m','MarkerSize',12)
axis([-0.3 0.3 0 1500]); set(gca,'YDir','reverse');

switch thd
case 0
subplot(2,3,4), hold on; grid on; box on; ylabel('PROFONDEUR (M)');
case 15
subplot(2,3,5), hold on; grid on; box on;
case 22
subplot(2,3,6), hold on; grid on; box on;
end

plot(TPP(2:NC-1),ZC(2:NC-1),'vr','MarkerSize',12), plot(TSS(2:NC-1),ZC(2:NC-1),'vb','MarkerSize',12)
plot(TPS(2:NC-1),ZC(2:NC-1),'vg','MarkerSize',12), plot(TSP(2:NC-1),ZC(2:NC-1),'vm','MarkerSize',12)
axis([-0.5 1.2 0 1500]); set(gca,'YDir','reverse');

% calcul pour frequences utilisees par ondelette
iom = 0;
for fo = df:df:(length(FG)-1)*df 
om = 2*pi*fo; iom = iom+1;

DM = eye(4); % matrice de transfert descendantes P et S - montantes P et S

EKPH = exp(i*om*qPH(2:NC)); CKPH = conj(EKPH); % propagation P dans couches
EKSH = exp(i*om*qSH(2:NC)); CKSH = conj(EKSH); % propagation S dans couches

% matrice transfert a travers couches et stratification

for ic = 3:NC-2
    
div = TPP(ic)*TSS(ic)-TPS(ic)*TSP(ic);
DP = [TSS(ic)*CKPH(ic) , -TSP(ic)*CKSH(ic) , (TSP(ic)*rPS(ic)-TSS(ic)*rPP(ic))*EKPH(ic) , (TSP(ic)*rSS(ic)-TSS(ic)*rSP(ic))*EKSH(ic)]/div;
DS = [TPS(ic)*CKPH(ic) , -TPP(ic)*CKSH(ic) , (TPP(ic)*rPS(ic)-TPS(ic)*rPP(ic))*EKPH(ic) , (TPP(ic)*rSS(ic)-TPS(ic)*rSP(ic))*EKSH(ic)]/(-div);
MP = RPP(ic)*DP + RSP(ic)*DS + [0 , 0 , tPP(ic)*EKPH(ic) , tSP(ic)*EKSH(ic)];
MS = RPS(ic)*DP + RSS(ic)*DS + [0 , 0 , tPS(ic)*EKPH(ic) , tSS(ic)*EKSH(ic)];

DM = DM*[DP ; DS ; MP ; MS];
end

DM = DM*[1 , 0 ; 0 , 1 ; RPP(NC-1) , RSP(NC-1) ; RPS(NC-1) , RSS(NC-1)]; % interface avec 1/2 espace inferieur

% cas sismique reflexion terrestre
% onde incidente P d'amplitude 1 au sommet stratification
DP1 = DM(1:2,:)\[1 ; 0]; %DP1 dans couche du bas 
RHPPS = DM(3:4,:)*DP1;  % ondes montantes P et S au sommet de la stratification
RHPP = RHPPS(1); RHPS = RHPPS(2); %coef. reflexion sur stratification elastique

AA(iom+1) =  RHPP*EKPH(2)*EKPH(2)*FG(iom+1); % ondes P sous fond de l'eau
BB(iom+1) =  RHPS*EKPH(2)*EKSH(2)*FG(iom+1); % ondes S sous fond de l'eau

% onde incidente S d'amplitude 1 au sommet stratification
DS1 = DM(1:2,:)\[0 ; 1]; %DP1 dans couche du bas 
RHSPS = DM(3:4,:)*DS1;  % ondes montantes P et S au sommet de la stratification
RHSP = RHSPS(1); RHSS = RHSPS(2); %coef. reflexion sur stratification elastique

CC(iom+1) =  RHSP*EKSH(2)*EKPH(2)*FG(iom+1); % ondes P sous fond de l'eau
DD(iom+1) =  RHSS*EKSH(2)*EKSH(2)*FG(iom+1); % ondes S sous fond de l'eau

% cas sismique reflexion marine
% matrice transfert a travers couche et stratification sous eau
DP = [CKPH(2) , 0 , -rPP(2)*EKPH(2) , -rSP(2)*EKSH(2)]/TPP(2);
MP = RPP(2)*DP + [0 , 0, tPP(2)*EKPH(2) , tSP(2)*EKSH(2)];
DME = [DP ; MP]*DM;

DE1 = DME(1,:)\1; %DP1 et DS1 pour onde incidente P d'amplitude 1 dans l'eau
RHEE = DME(2,:)*DE1; %coef. reflexion sur stratification elastique sous eau
RHEPS = DM(3:4,:)*DE1; RHEP = RHEPS(1); RHES = RHEPS(2); % ondes montantes P et S dans couche sous eau

FF(iom+1) =  RHEE*FG(iom+1); % onde P au fond de l'eau reflechie par stratification 
GG(iom+1) =  RHEP*EKPH(2)*FG(iom+1); % ondes P sous fond de l'eau reflechie par stratification
HH(iom+1) =  RHES*EKSH(2)*FG(iom+1); % ondes S sous fond de l'eau reflechie par stratification

end % frequences

ZE = zeros(nt/2+1-length(FF)); 

FP = [AA,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTPPB = max(abs(TF)); if MTPPA > 1e-6, TPPB = TF/MTPPA(ith); else TPPB = TF; end
FP = [BB,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTPSB = max(abs(TF)); if MTPSA > 1e-6, TPSB = TF/MTPSA(ith); else TPSB = TF; end
FP = [CC,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTSPB = max(abs(TF)); if MTSPA > 1e-6, TSPB = TF/MTSPA(ith); else TSPB = TF; end
FP = [DD,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTSSB = max(abs(TF)); if MTSSA > 1e-6, TSSB = TF/MTSSA(ith); else TSSB = TF; end
FP = [FF,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTEEB = max(abs(TF)); if MTEEA > 1e-6, TEEB = TF/MTEEA(ith); else TEEB = TF; end
FP = [GG,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTEPB = max(abs(TF)); if MTEPA > 1e-6, TEPB = TF/MTEPA(ith); else TEPB = TF; end
FP = [HH,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = real(ifft(conj(TF))); MTESB = max(abs(TF)); if MTESA > 1e-6, TESB = TF/MTESA(ith); else TESB = TF; end

PPAB = TPPB-TPPA(ith,:); MPPAB = max(abs(PPAB)); if MPPAB > 1e-6, PPAB = PPAB/MPPAB; end
PSAB = TPSB-TPSA(ith,:); MPSAB = max(abs(PSAB)); if MPSAB > 1e-6, PSAB = PSAB/MPSAB; end
SPAB = TSPB-TSPA(ith,:); MSPAB = max(abs(SPAB)); if MSPAB > 1e-6, SPAB = SPAB/MSPAB; end
SSAB = TSSB-TSSA(ith,:); MSSAB = max(abs(SSAB)); if MSSAB > 1e-6, SSAB = SSAB/MSSAB; end
EEAB = TEEB-TEEA(ith,:); MEEAB = max(abs(EEAB)); if MEEAB > 1e-6, EEAB = EEAB/MEEAB; end
EPAB = TEPB-TEPA(ith,:); MEPAB = max(abs(EPAB)); if MEPAB > 1e-6, EPAB = EPAB/MEPAB; end
ESAB = TESB-TESA(ith,:); MESAB = max(abs(ESAB)); if MESAB > 1e-6, ESAB = ESAB/MESAB; end

switch thd
    case 0, figure(3)
    case 15, figure(4)
    case 22, figure (5)
end


subplot(1,5,1)
for ic = 1:NC-3  plot([-1,-0.5],[TTPP(ic),TTPP(ic)],'r:'); plot([0.5,1],[TTPS(ic)-TTPS(1)+TTPP(1),TTPS(ic)-TTPS(1)+TTPP(1)],'g:'); end
for ic = 2:NC-3  plot([1,1.5],[TTPS(ic)-TTPS(2)+TTPP(2),TTPS(ic)-TTPS(2)+TTPP(2)],'g:'); end 
plot(TPPB,tt,'r:'), plot(2+PPAB,tt,'r')
axis([-1 3 0 1.4]); set(gca,'YDir','reverse');

subplot(1,5,2)
for ic = 1:NC-3  plot([-1,-0.5],[TTSS(ic),TTSS(ic)],'b:'); plot([0.5,1],[TTPS(ic)-TTPS(1)+TTSS(1),TTPS(ic)-TTPS(1)+TTSS(1)],'g:'); end
for ic = 2:NC-3  plot([1,1.5],[TTPS(ic)-TTPS(2)+TTSS(2),TTPS(ic)-TTPS(2)+TTSS(2)],'g:'); end
plot(TSSB,tt,'b:'), plot(2+SSAB,tt,'b') 
axis([-1 3 0 1.4]); set(gca,'YDir','reverse');

subplot(1,5,3)
plot(10*TEEB,tt,'k:'), plot(2+EEAB,tt,'k')
legend('TRACE A','TRACE B','TRACE B - TRACE A' ,['p = sin( ' int2str(thd) '° ) /1500 S/M']);
axis([-1 3 0 1.4]); set(gca,'YDir','reverse');

subplot(1,5,4)
for ic = 1:NC-3  plot([-1,-0.5],[TTPP(ic),TTPP(ic)],'r:'); plot([0.5,1],[TTPP(1)+TTPP(ic),TTPP(1)+TTPP(ic)],'r:'); 
plot([1,1.5],[TTPP(2)+TTPP(ic),TTPP(2)+TTPP(ic)],'r:'); plot([2.5,3],[TTPP(4)+TTPP(ic),TTPP(4)+TTPP(ic)],'r:'); end
plot(TEPB,tt,'r:'), plot(2+EPAB,tt,'r') 
axis([-1 3 0 1.4]); set(gca,'YDir','reverse'); 

subplot(1,5,5)
for ic = 1:NC-3  plot([-1,-0.5],[TTPS(ic),TTPS(ic)],'g:');  end
plot(TESB,tt,'g:'), plot(2+ESAB,tt,'g')
axis([-1 3 0 1.4]); set(gca,'YDir','reverse'); 

end % incidence