%HODOVZO dessine hodochrone et rayons quand un gradient de vitesse
% change de signe entre ZA<z<ZB, produisant une zone d'ombre 
% et un guide d'onde pour une source dans zone a faible vitesse
% voir : http://jmmeost.free.fr/Sismique/osis.html
% fig1  = http://jmmeost.free.fr/Sismique/hodovzo.gif  
% fig2  = http://jmmeost.free.fr/Sismique/hodovzg.gif  

clear all, clf

dz = 1; V0 = 1000; a = 1; 
ZA = 250; ZB = 500; XM = 5000; ZM = 2000; 
z = 0:dz:ZM; V = V0+a*z; VA = V0+a*ZA; 

% la source ponctuelle est a l'origine
figure(1)
subplot(2,1,1), hold on, grid on
title('HODOCHRONE POUR V = V0+AZ V0 = 1000 M/S A = +1 /S (Z < 250 M et Z > 500 M) A = -1 /S (250 < Z < 500 M)')
ylabel('TEMPS (MS)  VITESSE (M/S)')
set(gca,'Ydir','reverse'),axis([0,XM,0,XM])

subplot(2,1,2), hold on, grid on
title('ZONE D"OMBRE')
ylabel('PROFONDEUR (M)'), xlabel('DISTANCE (M)')
set(gca,'Ydir','reverse','DataAspectRatio',[1,1,1]),axis([0,XM,0,ZM-500])

% rayons avec Zmax<=ZA
for p = [1/V0:-0.00002:1/V(ZA),1/V(ZA)]
th0 = asin(p*V0); VZ = 1/p; ZZ = (VZ-V0)/a; 
R = VZ/a; XZ = R*cos(th0); TZ = -log(tan(th0/2))/a*1000;
% point et t courant du rayon entre 0<z<Zmax
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 = 0:.5:4 s
itf = 1;
for nt = 1:8
   itf = [itf, min(find((T-nt*500)>=0))];
end

subplot(2,1,1), plot(2*XZ,2*TZ,'b+',2*XZ,VZ,'co')
subplot(2,1,2), plot(X,Z,'b') 
plot(X(itf),Z(itf),'bo'), plot(XZ,ZZ,'m+')
end


% rayons traversant la tranche ZA<Z<ZB construits par mise bout a bout 
% d'arcs de cercle dont la courbure est vers le bas entre ZA et ZB 
for p = [1/V(ZA+1):-.00001:1/V(450),1/V(460):-.00002 :1/V(ZM)]
th0 = asin(p*V0); VZ = 1/p; ZZ = (VZ-V0)/a; 
R = VZ/a; XZ = R*cos(th0); TZ = -log(tan(th0/2))/a*1000;
thA = asin(p*VA); XA = R*(cos(th0)-cos(thA)); 
TA = log(tan(thA/2))/a*1000+TZ;
% point et t courant du rayon pour gradient constant
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;
Xa = X(find(Z)<=ZA); Za = Z(find(Z)<=ZA); Ta = T(find(Z)<=ZA);
Xb = X(find(Z)>ZA); Zb = Z(find(Z)>ZA); Tb = T(find(Z)>ZA);
% rayon entre 0<z<Zmax
X = [Xa,2*XA-fliplr(Xa),2*XA+Xa,Xb+2*XA]; 
Z = [Za,ZB-fliplr(Za),ZB+Za,Zb+ZB];
T = [Ta,2*TA-fliplr(Ta),2*TA+Ta,Tb+2*TA];
% tout le rayon
X = [X,2*XA+XZ,2*XZ+4*XA-fliplr(X)]; Z = [Z,ZB+ZZ,fliplr(Z)];
T = [T,TZ,2*TZ+4*TA-fliplr(T)];
% indice de position du front d'onde a t = .5:.5:4 s
itf = min(find((T-500)>=0));
for nt = 2:8
   itf = [itf, min(find((T-nt*500)>=0))];
end

subplot(2,1,1)
if p>1/V(450) % branche retrograde de l'hodochrone
plot(2*XZ+4*XA,2*TZ+4*TA,'r+',2*XZ+4*XA,VZ,'ro')
else % branche prograde de l'hodochrone
plot(2*XZ+4*XA,2*TZ+4*TA,'g+',2*XZ+4*XA,VZ,'go')
end
subplot(2,1,2)
if p>1/V(450) % rayons correspondant a branche retrograde
   plot(X,Z,'r'), plot(X(itf),Z(itf),'ro') 
else % rayons correspondant a branche prograde
   plot(X,Z,'g'), plot(X(itf),Z(itf),'go') 
end
plot(XZ+2*XA,ZZ+2*ZA,'m+');
end

subplot(2,1,1), plot(0,0,'b+',0,1000,'co')
legend('T(X(p))','1/p = V0/sinth0 = V(Zm)')

subplot(2,1,2), plot([0,XM],[ZA,ZA],'k',[0,XM],[ZB,ZB],'k')
legend('Rayons','Fronts Dt=.5s','Zm(p)')


% la source est dans la zone a basse vitesse en z=ZB
figure(2)
subplot(2,1,1), hold on, grid on
title('HODOCHRONE POUR V = V0+AZ V0 = 1000 M/S A = +1 /S (Z < 250 M et Z > 500 M) A = -1 /S (250 < Z < 500 M)')
ylabel('TEMPS (MS)  VITESSE (M/S)')
set(gca,'Ydir','reverse'),axis([0,XM,0,XM])

subplot(2,1,2), hold on, grid on
title('ZONE D"OMBRE ET GUIDE D"ONDE')
ylabel('PROFONDEUR (M)'), xlabel('DISTANCE (M)')
set(gca,'Ydir','reverse','DataAspectRatio',[1,1,1]),axis([0,XM,0,ZM-500])

% rayons avec Zmax dans zone basse vitesse
for p = [1/V0-.00002:-0.00002:1/V(ZA),1/V(ZA)]
th0 = asin(p*V0); VZ = 1/p; ZZ = (VZ-V0)/a+ZB; 
R = VZ/a; XZ = R*cos(th0); TZ = -log(tan(th0/2))/a*1000;
% point et t courant du rayon partant vers le bas
pv = p*V; pv = pv(find(pv<1)); th = asin(pv); 
X = [XZ-R*cos(th)]; Z = [z(find(pv<1))]+ZB;
T = log(tan(th/2))/a*1000+TZ;
% tout le rayon partant vers le bas
X = [X,XZ,2*XZ-fliplr(X)]; Z = [Z,ZZ,fliplr(Z)]; 
T = [T,TZ,2*TZ-fliplr(T)];
% rayon partant vers le bas et le haut
X = [X,X]; Z = [Z,2*ZB-Z]; T = [T,T];
% trajets successifs dans guide d'onde
for tra = 2:2:9 
   X = [X,X+tra*XZ]; Z = [Z,Z]; T = [T,T+tra*TZ]; 
end
% indice de position du front d'onde a t = .5:.5:2 s
itf = min(find((T-500)>=0));
for nt = 2:4
   itf = [itf, min(find((T-nt*500)>=0))];
end 

subplot(2,1,1), %plot(2*XZ,2*TZ,'r+',2*XZ,VZ,'ro')
subplot(2,1,2), plot(X,Z,'r'),plot(XZ,2*ZB-ZZ,'m+')
plot(X(itf),Z(itf),'ro',X(itf),2*ZB-Z(itf),'ro'),plot(XZ,ZZ,'m+')

end

% rayons sortant de zone a basse vitesse 
for p = [1/V(ZA)-.00002:-.00002:1/V(ZM),.00032:-.00002:.00002]
th0 = asin(p*V0); VZ = 1/p; ZZ = ZB-(VZ-V0)/a; 
R = VZ/a; XZ = R*cos(th0); TZ = -log(tan(th0/2))/a*1000;
thA = asin(p*VA); XA = R*(cos(th0)-cos(thA)); 
TA = log(tan(thA/2))/a*1000+TZ;
% point et t courant du rayon 
pv = p*V; pv = pv(find(pv<1)); th = asin(pv); 
X = [XZ-R*cos(th)]; T = log(tan(th/2))/a*1000+TZ;

% rayon partant vers le haut
Z = ZB-[z(find(pv<1))];
Xa = X(find(Z)<=ZA); Za = Z(find(Z)<=ZA); Ta = T(find(Z)<=ZA);
% rayon entre 0<z<Zmax
XH = [Xa,2*XA-fliplr(Xa)];  ZH = [Za,ZB-fliplr(Za)]; 
TH = [Ta,2*TA-fliplr(Ta)]; 
% indice de position du front d'onde a t = .5 s
itfh = [min(find((TH-500)>=0)),min(find((TH-1000)>=0))];    

% rayons  partant vers le bas
Z = 2*ZB-Z;
% le rayon pour z>ZB
XD = [X,XZ,2*XZ-fliplr(X)]; ZD = [Z,2*ZB-ZZ,fliplr(Z)]; 
TD = [T,TZ,2*TZ-fliplr(T)];
% rayon pour z>ZB et z<ZB
XD = [XD,2*XZ+XH]; ZD = [ZD,ZH]; TD = [TD,2*TZ+TH];
% indice de position du front d'onde a t = 0.5:.5:4 s
itfd = min(find((TD-500)>=0));
for nt = 2:8
   itfd = [itfd, min(find((TD-nt*500)>=0))];
end 

subplot(2,1,1)
plot(2*XA,2*TA,'b+',2*XA,VZ,'co')
plot(XD(length(XD)),TD(length(TD)),'g+',XD(length(XD)),VZ,'go')
subplot(2,1,2)
plot(XH,ZH,'b'), plot(XH(itfh),ZH(itfh),'bo')
plot(XD,ZD,'g'), plot(XD(itfd),ZD(itfd),'go'),  plot(XZ,2*ZB-ZZ,'m+')
end

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

subplot(2,1,2), legend('Rayons','Zm(p)','Fronts Dt=.5s')
plot([0,XM],[ZA,ZA],'k',[0,XM],[ZB,ZB],'k')