%HODOSTRAT dessine hodochrone pour ondes reflechies et coniques en milieu stratifié
%
% Documenter l'effet de variations de vitesse et épaisseurs de couches sur 
% l'hodochrone sur les réflexions et les ondes coniques (existence,
%  observation en 1ere arrivée)

clear all
cl = ['b' 'r' 'g' 'm' 'y']; % couleur des interfaces

%VALEURS NUMERIQUES*******************************************************
xmax = 10000; % longueur dispositif enregistrement en m
dg = 10; xg = 0:dg:xmax; % intervalle et abscisse récepteurs en m
tmax = 5; % temps max d'enregistrement en s

dp = .00001; % intervalle paramètre pour tracé de rayons en s/m

% modele vitesses à 5 couches
HC =  [ 250  500  750 1000  500]'; % epaisseur couches en m
VP =  [1500 2000 3000 4000 5000]'; % vitesses P en m/s
%*************************************************************************

nc = length(HC);  ZC = cumsum(HC); zmax = max(ZC); % nombre couches
tip = 2*HC./VP; % temps double vertical P dans chaque couche

figure (1), clf

% MODELE VITESSE**********************************************************
subplot(2,2,1); hold on; grid on; box on;
stairs([VP(1) ; VP],[0 ; ZC],'k','Linewidth',4); % V(z)
for ic = 1:nc-1, plot([0 xmax],ZC(ic)+[0 0],cl(ic),'Linewidth',4); end % interfaces
xlabel('VITESSES (M/S)');  ylabel('PROFONDEUR (M)'); title('MODELE VITESSE'); 
axis([0 xmax 0 zmax]); set(gca,'YDir','reverse','Color',[1. 0.8 0.5])


% HODOCHRONE REFLEXIONS***************************************************
% temps arrivée onde directe et reflechie dans couche 1
subplot(2,2,2), hold on, box on, grid on, 
plot([0 xmax],[0 xmax/VP(1)],'b','LineWidth',3) % onde directe
plot(xg,sqrt(tip(1)^2+(xg/VP(1)).^2),'b','LineWidth',3) % onde réfléchie sur interface 12
xlabel('X (M)'), ylabel('T (S)'), title('TEMPS REFLECHIS ET CONIQUES')
axis([0 xmax 0 tmax]), set(gca,'YDir','reverse','Color',[1. 0.8 0.5])

% ondes P reflechies aux interfaces
pc = [0:dp:1/VP(1) 1/VP(1)]; % parametres pour incidence entre 0 et 90° dans couche 1
th1 = asin(pc*VP(1)); % angles incidence dans couche 1
xp = 2*HC(1)*tan(th1); % distances horizontale dans couche 1
tp = tip(1)./cos(th1); % temps trajet dans couche 1
for ic = 2:nc-1 % boucle sur interfaces
pc = [0:dp:1/max(VP(1:ic)) 1/max(VP(1:ic))]; npc = length(pc); % parametres pour incidence dans couche i
thi = asin(pc*VP(ic)); % angle incidence dans couche i
xp = xp(1:npc)+2*HC(ic)*tan(thi); % distance horizontale cumulée dans couche 1 à i
tp = tp(1:npc)+tip(ic)./cos(thi); % temps trajet cumulé dans couche 1 à i
plot(xp,tp,cl(ic),'LineWidth',3)
end

% HODOCHRONE ET RAYONS CONIQUES*******************************************
% interfaces donnant ondes coniques
icc = 0;
for ic = (find((VP-circshift(VP,1))>0))'-1 % interfaces avec V croissante
if VP(ic+1) > max(VP(1:ic)) % et V > V dans toutes les couches au-dessus
if icc == 0, icc = ic; else icc = [icc ic]; end % indices des interfaces donnant ondes P conique
end
end
ncc = length(icc); % nombre interfaces donnant coniques

subplot(2,2,4), hold on, box on, grid on, 
xlabel('X (M)');  ylabel('PROFONDEUR (M)'); title('RAYONS CONIQUES'); 
axis([0 xmax 0 zmax]); set(gca,'YDir','reverse','DataAspectRatio',[1 1 1],'Color',[1. 0.8 0.5])

% ondes coniques
for ic = icc % boucle sur interfaces donnant ondes coniques
thc = asin(VP(1:ic)/VP(ic+1)); % loi de Snell
xc = 2*sum(HC(1:ic).*tan(thc)); % distance critique
tc = sum(tip(1:ic)./cos(thc)); % temps arrivée à xc
tm = tc+(xmax-xc)/VP(ic+1); % temps arrivée à xmax

subplot(2,2,2), % temps arrivée
plot([xc xmax],[tc tm],[cl(ic+1) '--'],'MarkerSize',8,'LineWidth',3)

subplot(2,2,4) % rayons
if ic == 1,  % conique sous interface 1
plot([0 xc/2 xc],[0 ZC(ic) 0],cl(ic+1),xmax-[xc/2 0],[ZC(ic) 0],cl(ic+1),[xc/2 xmax-xc/2],ZC(ic)+[0 0],cl(ic+1),'LineWidth',3), end
if ic > 1 % conique sous interface i
xi = HC(1:ic-1).*tan(thc(1:ic-1)); xii = cumsum(xi); % intersection rayon - interface 1 2 et 3
plot([0 xii' xc/2 xc-xii(ic-1:-1:1)' xc],[0 ZC(1:ic)' ZC(ic-1:-1:1)' 0],cl(ic+1),xmax-[xc/2 xii(ic-1:-1:1)' 0],[ZC(ic:-1:1)' 0],cl(ic+1),[xc/2 xmax-xc/2],ZC(ic)+[0 0],cl(ic+1),'LineWidth',3)
end
end