%DIRECO dessine hodochrones et fronts d'onde direct, reflechi, conique 
% et transmis pour source ponctuelle a distance d d'une interface plane
% voir : http://jmmeost.free.fr/Sismique/osis.html
% fig1 = http://jmmeost.free.fr/Sismique/hodo.gif  
% fig2 = http://jmmeost.free.fr/Sismique/hodop.gif

clear all

V1 = 1000; V2 = 2000; d = 500; alpha = 15;
XM = 4000; TM = 3000; x = -XM:1:XM; 

figure(1), clf
subplot(2,1,2), hold on, grid on
title(['HODOCHRONE INTERFACE HORIZONTALE V1 = ' num2str(V1) ' M/S , V2 = ' num2str(V2) ' M/S , D = ' num2str(d) ' M'])
xlabel('DISTANCE (M)'), ylabel('TEMPS (MS)')

% theta , x , t pour incidence critique
thc = asin(V1/V2); xc = 2*d*tan(thc); tc = 2*d/cos(thc)/V1*1000;

% t direct en ms
td = abs(x)/V1*1000;
plot(x,td,'b','Linewidth',3)

% t reflechi en ms
tr = sqrt(4*d^2+x.^2)/V1*1000;
plot(x,tr,'r','Linewidth',3)

% t conique droite et gauche en ms
tcd = (x(find(x>xc))/V2+2*d*cos(thc)/V1)*1000;
tcg = (-x(find(x<-xc))/V2+2*d*cos(thc)/V1)*1000;
plot(x(find(x>xc)),tcd,'g',x(find(x<-xc)),tcg,'g','Linewidth',3)

plot(xc,TM,'g+',-xc,TM,'g+',-XM,tc,'g+',xc,tc,'g+',-xc,tc,'g+','Markersize',10)
text(xc,TM-150,'Xc','Fontsize',12), text(-xc,TM-150,'Xc','Fontsize',12), text(-XM+100,tc,'Tc','Fontsize',12)

legend('DIRECTE','REFLECHIE','CONIQUE'), 

set(gca,'YDir','reverse','DataAspectRatio',[1,1,1],'Color',[1. 0.8 0.5]), axis([-XM,XM,0,TM])


subplot(2,1,1), hold on, , grid on
title('FRONTS ONDES DIRECT REFLECHI CONIQUE TRANSMIS A T = 0.75 , 1 , 1.5 , 2 , 2.5 S ')
xlabel('DISTANCE (M)'), ylabel('PROFONDEUR (M)')
set(gca,'YDir','reverse','DataAspectRatio',[1,1,1],'Color',[1. 0.8 0.5]), axis([-XM,XM,0,TM])
plot(xc,0,'g+',-xc,0,'g+','Markersize',10)
plot([-XM,XM],[d,d],'k')

for Dt = [0.75,1.:.5:2.5] ;
% front d'onde direct spherique
z = V1^2*Dt^2-x.^2; 
zs = sqrt(z(find(z>=0))); xs = x(find(z>=0));
zfd = zs(find(zs<=d)); xfd = xs(find(zs<=d));
plot(xfd(find(xfd>=0)),zfd(find(xfd>=0)),'b','Linewidth',3)
plot(xfd(find(xfd<=0)),zfd(find(xfd<=0)),'b','Linewidth',3)

% front d'onde reflechi issu source virtuelle
zfr = 2*d-zs(find((zs>=d)&(zs<=2*d))); 
xfr = xs(find((zs>=d)&(zs<=2*d)));
plot(xfr(find(xfr>=0)),zfr(find(xfr>=0)),'r','Linewidth',3)
plot(xfr(find(xfr<=0)),zfr(find(xfr<=0)),'r','Linewidth',3)

% front d'onde conique entre front reflechi et interface 
xfrc = Dt*V1*sin(thc); zfrc = 2*d-Dt*V1*cos(thc); 
xfc = (Dt-d*cos(thc)/V1)*V2;
plot([xfrc,xfc],[zfrc,d],'g','Linewidth',3)
plot([-xfrc,-xfc],[zfrc,d],'g','Linewidth',3)

% front d'onde transmis par trace rayons
th1 = 0:0.01:thc; th2 = asin(V2/V1*sin(th1));
xft = d*tan(th1)+(Dt-d./cos(th1)/V1)*V2.*sin(th2); xft = [xft,xfc];
zft = d+(Dt-d./cos(th1)/V1)*V2.*cos(th2); zft = [zft,d];
plot(xft,zft,'m','Linewidth',3)
plot(-xft,zft,'m','Linewidth',3)

end


figure(2), clf
subplot(2,1,2), hold on, grid on
title(['HODOCHRONE INTERFACE PENTEE ALPHA = ' num2str(alpha) ' ° ,  V1 = ' num2str(V1) ' M/S , V2 = ' num2str(V2) ' M/S , D = ' num2str(d) ' M'])
ylabel('TEMPS (MS)')

% interface pendage alpha intersecte surface en xmax
al = alpha*pi/180; xmax = d/sin(al); tmax = xmax/V1*1000;

% t direct
td = [abs(x(find(x<xmax)))/V1*1000,(xmax/V1+(x(find(x>=xmax))-xmax)/V2)*1000];
plot(x,td,'b','Linewidth',3)

% t reflechi issu source virtuelle, minimum en trm 
xv = 2*d*sin(al); zv = 2*d*cos(al);
tr = (sqrt(zv^2+(x-xv).^2)/V1)*1000; trm = zv/V1*1000;
plot(x(find(x<xmax)),tr(find(x<xmax)),'r','Linewidth',3)

% x, t critique aval, amont
xcv = 2*d*(cos(al)*tan(thc+al)-sin(al)); 
tcv = (xcv*sin(thc+al)+2*d*cos(thc))/V1*1000;
xcm = 2*d*(cos(al)*tan(thc-al)+sin(al)); 
tcm = (xcm*sin(thc-al)+2*d*cos(thc))/V1*1000;
% t conique aval , amont
tcav = (-x(find(-x>xcv))*sin(thc+al)+2*d*cos(thc))/V1*1000;
tcam = (x(find((x>xcm)&(x<xmax)))*sin(thc-al)+2*d*cos(thc))/V1*1000;
% t transmis en x > xmax 
tt = ((xmax*sin(thc-al)+2*d*cos(thc))/V1+(x(find(x>=xmax))-xmax)/V2)*1000;
plot(x(find((x>xcm)&(x<xmax))),tcam,'g','Linewidth',3)
plot(x(find(x>=xmax)),tt,'m','Linewidth',3)
plot(x(find(-x>xcv)),tcav,'g','Linewidth',3)

legend('DIRECTE','REFLECHIE','CONIQUE','TRANSMISE')
plot(-xcv,TM,'g+',-xcv,tcv,'g+',-XM,tcv,'g+','Markersize',10), 
text(-xcv,TM-150,'Xc aval','Fontsize',12), text(-XM+100,tcv,'Tc aval','Fontsize',12)
plot(xcm,TM,'g+',xcm,tcm,'g+',-XM,tcm,'g+','Markersize',10), 
text(xcm,TM-150,'Xc amont','Fontsize',12), text(-XM+100,tcm,'Tc amont','Fontsize',12)
plot(xv,trm,'r+',xv,TM,'r+','Markersize',10)
text(xv,TM-150,'Xv','Fontsize',12)
plot(xmax,tmax,'k+','Markersize',10)

set(gca,'YDir','reverse','DataAspectRatio',[1,1,1],'Color',[1. 0.8 0.5]), axis([-XM,XM,0,TM])


subplot(2,1,1), hold on, grid on
title('FRONTS ONDES DIRECT REFLECHI CONIQUE TRANSMIS A T = 0.75 , 1 , 1.5 , 2 , 2.5 S ')
xlabel('DISTANCE (M)'), ylabel('PROFONDEUR (M)')
set(gca,'YDir','reverse','DataAspectRatio',[1,1,1],'Color',[1. 0.8 0.5]), axis([-XM,XM,0,TM])

for Dt = [0.75,1.:.5:2.5] ;
% front d'onde direct   
z = V1^2*Dt^2-x.^2; 
zs = sqrt(z(find(z>=0))); xs = x(find(z>=0));
zfd = zs(find(zs<=(xmax-xs)*tan(al))); 
xfd = xs(find(zs<=(xmax-xs)*tan(al)));
plot(xfd(find(xfd>=0)),zfd(find(xfd>=0)),'b','Linewidth',3)
plot(xfd(find(xfd<=0)),zfd(find(xfd<=0)),'b','Linewidth',3)

% front d'onde reflechi issu source virtuelle
zr = V1^2*Dt^2-(x-xv).^2; 
zs = zv-sqrt(zr(find(zr>=0))); xs = x(find(zr>=0));
zfr = zs(find(zs<=(xmax-xs)*tan(al))); 
xfr = xs(find(zs<=(xmax-xs)*tan(al)));
plot(xfr(find(xfr>=0)),zfr(find(xfr>=0)),'r','Linewidth',3)
plot(xfr(find(xfr<=0)),zfr(find(xfr<=0)),'r','Linewidth',3)

% front d'onde conique aval et amont entre front reflechi et interface 
l1 = d/cos(thc); l2 = (Dt-l1/V1)*V2;
zfrcav = l1*cos(thc-al)-(V1*Dt-l1)*cos(thc+al); 
xfrcav = l1*sin(thc-al)+(V1*Dt-l1)*sin(thc+al);
zfcav = l1*cos(thc-al)+l2*sin(al); 
xfcav = l1*sin(thc-al)+l2*cos(al);
plot([-xfrcav,-xfcav],[zfrcav,zfcav],'g','Linewidth',3)

zfrcam = l1*cos(thc+al)-(V1*Dt-l1)*cos(thc-al); 
xfrcam = l1*sin(thc+al)+(V1*Dt-l1)*sin(thc-al);
zfcam = l1*cos(thc+al)-l2*sin(al); 
xfcam = l1*sin(thc+al)+l2*cos(al);
plot([xfrcam,xfcam],[zfrcam,zfcam],'g','Linewidth',3)

% front d'onde transmis
th1 = 0:0.01:thc; th2 = asin(V2/V1*sin(th1));
xftam = d./cos(th1).*sin(th1+al)+(Dt-d./cos(th1)/V1)*V2.*sin(th2+al); 
xftam = [xftam,xfcam];
zftam = d./cos(th1).*cos(th1+al)+(Dt-d./cos(th1)/V1)*V2.*cos(th2+al); 
zftam = [zftam,zfcam];
plot(xftam,zftam,'m','Linewidth',3)

xftav = d./cos(th1).*sin(th1-al)+(Dt-d./cos(th1)/V1)*V2.*sin(th2-al); 
xftav = [xftav,xfcav];
zftav = d./cos(th1).*cos(th1-al)+(Dt-d./cos(th1)/V1)*V2.*cos(th2-al); 
zftav = [zftav,zfcav];
plot(-xftav,zftav,'m','Linewidth',3)

end

plot(xcm,0,'g+',-xcv,0,'g+'), plot(0,0,'b+',xv,zv,'r+','Markersize',10)
plot([xmax,-XM],[0,(XM+xmax)*tan(al)],'k')