%HODOVZ dessine l'hodochrone, les rayons et fronts d'onde 
% dans un milieu a gradient de vitesse constant croissant
% pour une source en z = 0 et en z = ZS
% voir : http://jmmeost.free.fr/Sismique/osis.html
% fig  = http://jmmeost.free.fr/Sismique/hodovz.gif  

clear all, clf

dz = 1; V0 = 1000; a = 1;
XM = 5000; TM = 3000; ZM = 3000; ZS = 1000;

z = 0:dz:ZM; V = V0+a*z; VS = V0+a*ZS;

% source ponctuelle a l'origine

subplot(2,2,1), hold on, grid on
title('HODOCHRONE ET VITESSE APPARENTE POUR V = 1000+Z M/S')
ylabel('TEMPS (MS) , VITESSE (M/S)')
set(gca,'Ydir','reverse'),axis([0,XM,0,TM])

% rayon d'incidence theta0(t) , hodochrone et vitesses apparente
t = .01:.1:3.1;
th0 = 2*atan(exp(-a*t/2)); p = sin(th0)/V0; VZ = 1./p;
R = VZ/a; X = 2*R.*cos(th0); 
plot(X,t*1000,'b',X,VZ,'r')

subplot(2,2,3), hold on, grid on
title('FRONTS D"ONDES A T = (.5n) S ET RAYONS')
ylabel('PROFONDEUR (M)'), xlabel('DISTANCE (M)')
set(gca,'Ydir','reverse','DataAspectRatio',[1,1,1]),axis([0,XM,0,ZM])

for t = .5:.5:5
% rayon emergeant a t sur hodochrone
th0 = 2*atan(exp(-a*t/2)); p = sin(th0)/V0; VZ = 1/p; ZZ = (VZ-V0)/a; 
R = VZ/a; XZ = R*cos(th0); TZ = t/2*1000;
subplot(2,2,1), plot(0,0,'b+',2*XZ,2*TZ,'b+',2*XZ,VZ,'ro')

% point et t courant du rayon a la descente
pv = p*V; pv = pv(find(pv<1)); th = asin(pv);
X = [XZ-R*cos(th)]; Z = [z(find(pv<1))]; T = log(tan(th/2))/a*1000+TZ;
% tout le rayon
X = [X,XZ,2*XZ-fliplr(X)]; Z = [Z,ZZ,fliplr(Z)]; 
T = [T,TZ,2*TZ-fliplr(T)]; 
% indice de position du front d'onde a t = .5:.5:3
itf = 1;
for nt = 1:6
   itf = [itf, min(find((T-nt*500)>=0))];
end

%dessine le rayon et la position du front d'onde a t sur rayon
subplot(2,2,3), plot(X,Z,'b'), plot(X(itf),Z(itf),'ro')
    
% front d'onde a t est le cercle de centre 0,ZT et rayon XT
tht = 2*atan(exp(-a*t)); 
RT = V0/sin(tht)/a; XT = RT*cos(tht); ZT = RT-V0/a;
XF = XT^2-(z-ZT).^2; ZF = z(find(XF>=0)); XF = sqrt(XF(find(XF>=0)));

plot(XF,ZF,'r'), plot(0,ZZ,'m+',XZ,ZZ,'m+')

end

subplot(2,2,1)
legend('T(X(p))','1/p = V0/sinth0 = V(Zm) ')

subplot(2,2,3)
legend('Rayons','Fronts','Fronts','Zmax')

% source en x = 0, z = ZS

subplot(2,2,2), hold on, grid on
title('SOURCE EN ZS = 1000M')
set(gca,'Ydir','reverse'),axis([0,XM,0,TM])

t = .001:.1:3;
ths = 2*atan(exp(-a*t/2)); p = sin(ths)/VS; VZ = 1./p;
R = VZ/a; X = 2*R.*cos(ths);
th0 = asin(p*V0); X0 = R.*(cos(th0)-cos(ths)); 
T0 = log(tan(ths/2)./tan(th0/2))/a*1000;
plot(X0,T0,'b',X0,VZ,'r')
plot(X0+X,T0+t*1000,'b',X0+X,VZ,'r')

subplot(2,2,4), hold on, grid on
title('FRONTS D"ONDES A T = (.2, .5+.5n) S ET RAYONS')
xlabel('DISTANCE (M)')
set(gca,'Ydir','reverse','DataAspectRatio',[1,1,1]),axis([0,XM,0,ZM])

for t = [.2,.5:.5:3,4,5.5,10]
% rayon passant a t en ZS
ths = 2*atan(exp(-a*t/2)); p = sin(ths)/VS; VZ = 1/p; ZZ = (VZ-VS)/a; 
R = VZ/a; XZ = R*cos(ths); TZ = t/2*1000;

% point et t courant du rayon a la descente depuis ZS
pv = p*V; pv = pv(find(pv<1)); th = asin(pv);
X = [XZ-R*cos(th)]; Z = [z(find(pv<1))]; T = log(tan(th/2))/a*1000+TZ;
% rayon partant vers le haut par symetrie
XH = -fliplr(X); ZH = fliplr(Z); TH = -fliplr(T);
% tout le rayon partant vers le bas
X = [X,XZ,2*XZ+XH]; Z = [Z,ZZ+ZS,ZH]; T = [T,TZ,2*TZ+TH];
% indice de position du front d'onde a t 
itf = min(find((T-200)>=0)); itfh = min(find((TH-200)>=0));
for nt = 1:6
   itf = [itf, min(find((T-nt*500)>=0))];
   itfh = [itfh, min(find((TH-nt*500)>=0))];
end
% dessine hodochrone et vitesse apparente
subplot(2,2,2)
plot(XH(find(ZH<1)),TH(find(ZH<1)),'b+'), plot(XH(find(ZH<1)),VZ,'ro')
plot(X(find(Z<1)),T(find(Z<1)),'b+'), plot(X(find(Z<1)),VZ,'ro')
% dessine rayon
subplot(2,2,4) 
plot(X,Z,'b'), plot(X(itf),Z(itf),'ro')
plot(XH,ZH,'b'), plot(XH(itfh),ZH(itfh),'ro')
plot(XZ,ZS+ZZ,'m+')
 
% front d'onde a t est le cercle de centre 0,ZT et rayon XT
tht = 2*atan(exp(-a*t)); 
RT = VS/sin(tht)/a; XT = RT*cos(tht); ZT = ZS+RT-VS/a;
XF = XT^2-(z-ZT).^2; ZF = z(find(XF>=0)); XF = sqrt(XF(find(XF>=0)));
% dessine front d'onde
subplot(2,2,4), plot(XF,ZF,'r')

end