%STRATPSV calcule ondes harmoniques montante + descendante P et SV dans stratification  
% et trace sismique par TF inverse
% voir : http://jmmeost.free.fr/Sismique/osis.html#CS
% fig1  = http://jmmeost.free.fr/Sismique/mernord.gif
% fig2  = http://jmmeost.free.fr/Sismique/reflmernord.gif
% fig3  = http://jmmeost.free.fr/Sismique/stratifpsv0.gif
% fig4  = http://jmmeost.free.fr/Sismique/stratifpsv1.gif
% fig5  = http://jmmeost.free.fr/Sismique/stratifpsv3.gif

clear all

% modele vitesses et densite type Mer du Nord
HC =  [   0, 750, 350, 300, 100, 300,  20,  20,  20,  20,  20,  20,  20,  60, 200,  15,  45,  30,  15,  30, 300, 365];
VP =  [1500,1500,1860,2066,2049,2176,1950,2440,1950,2440,1950,2440,1950,2286,2751,2900,2770,2913,2900,2822,2900,3000];
VS =  [   0,   0, 950, 950, 950,1000, 950,1108, 950,1108, 950,1108, 950,1152,1470,1550,1614,1584,1550,1646,1550,1650];
RHO = [1000,1000,1600,1750,1900,2089,1800,2100,1800,2100,1800,2100,1800,2144,2219,2400,2200,2281,2400,2180,2400,2500];

NC = length(HC); ZC = cumsum(HC); % NC-1 couches de 2 (eau) a NC (1/2 espace inf) ; NC-1 interfaces de 1 (surface libre) a NC-1 
VPVS = VP./VS;  
IP = RHO.*VP; IP12 = 2*filter(ones(1,2)/2,1,IP); % impedances P, IP12 = IP1+IP2...
RP0 = -diff(IP)./IP12(2:NC); RP0(1)=1; TP0 = 2*IP(2:NC-1)./IP12(3:NC); % coef. refl. trans. P incidence verticale
IS = RHO.*VS; IS12 = 2*filter(ones(1,2)/2,1,IS); % impedances S, IS12 = IS1+IS2...
RS0 = -diff(IS)./IS12(2:NC); RS0(2)=1; TS0 = 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(VP,ZC,'r'), stairs(VS,ZC,'b'),  stairs(1000*(VPVS-1.7),ZC,'m'), stairs(RHO,ZC,'g'), stairs(VP,ZC,'r'); 
legend('VP','VS','(VP/VS-1.7)*1000','RHO'); xlabel('VITESSES (M/S) , DENSITE (KG/M3)');  
ylabel('PROFONDEUR (M)'); title('STRATIFICATION TYPE MER DU NORD'); 
axis([0 3500 0 3000]); set(gca,'YDir','reverse','Color',[1. 0.8 0.5]);

subplot(1,2,2); hold on; grid on; box on;
stairs(IP*1e-7,ZC,'r'); stairs(IS*1e-7,ZC,'b'); 
plot(RP0,ZC(1:NC-1),'^r','MarkerFaceColor','r'); plot(RS0(3:NC-1),ZC(3:NC-1),'^b','MarkerFaceColor','b');
plot(TP0,ZC(2:NC-1),'vr','MarkerFaceColor','r'); plot(TS0,ZC(2:NC-1),'vb','MarkerFaceColor','b');
legend('IP*1e-7','IS*1e-7','RPP(theta = 0)','RSS(eta = 0)','TPP(theta = 0)','TSS(eta = 0)'); 
xlabel('COEF. REFLEXION , IMPEDANCES (KG/M2S) , COEF. TRANSMISSION'); 
axis([-.4 1.2 0 3000]); set(gca,'YDir','reverse','Color',[1. 0.8 0.5]);

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

% angle incidence dans l'eau
for thd = 0:15:30
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 30
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 legend('RPP','RSS','RPS','RSP'); end
axis([-0.3 0.3 0 3000]); set(gca,'YDir','reverse','Color',[1. 0.8 0.5]); 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 30
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 legend('TPP','TSS','TPS','TSP'); end
axis([-0.5 1.2 0 3000]); set(gca,'YDir','reverse','Color',[1. 0.8 0.5]); xlabel('COEF. TRANSMISSION');


% 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

iom = 0;
for fo = df:df:(length(FG)-1)*df %frequences utilisees par ondelette
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 = ifft(conj(TF)); MTPP = max(abs(TF)); TPP = TF/MTPP;
FP = [BB,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = ifft(conj(TF)); MTPS = max(abs(TF)); TPS = TF/MTPS;
FP = [CC,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = ifft(conj(TF)); MTSP = max(abs(TF)); TSP = TF/MTSP;
FP = [DD,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = ifft(conj(TF)); MTSS = max(abs(TF)); TSS = TF/MTSS;
FP = [FF,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = ifft(conj(TF)); MTEE = max(abs(TF)); TEE = TF/MTEE;
FP = [GG,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = ifft(conj(TF)); MTEP = max(abs(TF)); TEP = TF/MTEP;
FP = [HH,ZE(1,:)]; FN = conj(fliplr(FP)); TF = [FP,FN(2:nt/2)]; TF = ifft(conj(TF)); MTES = max(abs(TF)); TES = TF/MTES;

switch thd
    case 0, figure(3)
    case 15, figure(4)
    case 30, 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'); end
plot(real(TPP),tt,'r'); 
title(['P INC. SANS EAU']);
axis([-1 1 0 nt*dt]); set(gca,'YDir','reverse'); xlabel('P REFL. SANS EAU'); legend('TTPP'); 
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'); end
plot(real(TSS),tt,'b'); 
title(['SV INC. SANS EAU']);
axis([-1 1 0 nt*dt]); set(gca,'YDir','reverse'); xlabel('SV REFL. SANS EAU'); legend('TTSS');

subplot(1,5,3), plot(0,-1,'.w'), hold on; grid on; box on;
plot(real(TEE),tt,'k');
title(['P INC. DANS EAU']); xlabel('P REFL. DANS EAU'); legend(['REFLEXION P-SV SUR STRATIFICATION TYPE MER DU NORD p = sin(' int2str(thd) '°)/1500 S/M']);
axis([-1 1 0 nt*dt]); set(gca,'YDir','reverse');

subplot(1,5,4), plot(0,-1,'.w'), hold on; grid on; box on;
for ic = 1:NC-3  plot([-1,-0.5],[TTPP(ic),TTPP(ic)],'r'); end
plot(real(TEP),tt,'r'); 
title(['P INC. DANS EAU']); xlabel('P REFL. SOUS EAU'); legend(['MAX = ' num2str(100*MTEP,'%3.2g')],'TTPP');
axis([-1 1 0 nt*dt]); set(gca,'YDir','reverse'); 

subplot(1,5,5), plot(0,-1,'.w'), hold on; grid on; box on;
for ic = 1:NC-3  plot([-1,-0.5],[TTPS(ic),TTPS(ic)],'g');  end
plot(real(TES),tt,'g');
title(['P INC. DANS EAU']); xlabel('SV REFL. SOUS EAU'); legend(['MAX = ' num2str(100*MTES,'%3.2g')],'TTPS');
axis([-1 1 0 nt*dt]); set(gca,'YDir','reverse'); 

end % incidence