%ISIS illustre les differences entre structure en profondeur et section sismique en temps
% voir : http://jmmeost.free.fr/Sismique/isis.html
% fig1 = http://jmmeost.free.fr/Sismique/plan.gif  
% fig3 = http://jmmeost.free.fr/Sismique/courbe.gif
% fig4 = http://jmmeost.free.fr/Sismique/concave.gif

clear all

V = 2000;
z0 = 500; t0 = 2*z0/V;
lx = 1000; lz = 1000;
x0 = 300; x1 = 700; dx = 100; 

% reflecteur plan
figure(1), clf

subplot(2,2,1), hold on, grid on
plot(0,0,'ko',0,0,'kx'), plot(0,0,'r','LineWidth',2), plot(0,0,'r',0,0,'g')
legend('Sources','Récepteurs','Réflecteur','Rayons réfléchis','Rayons diffractés',['V = ' num2str(V) ' M/S'])
% rayons diffractes
X = 0:dx:1000; NX = length(X);
for ix = 1:NX
    plot([X(ix) x0],[0,z0],'g',[X(ix) x1],[0,z0],'g')
end
% rayons reflechis
X = x0:dx:x1; NX = length(X);
for ix = 1:NX
    plot(X(ix),0,'ro',X(ix),0,'rx',[X(ix) X(ix)],[0,z0],'r')
end
% pas reflexion
X0 = [0:dx:x0-dx , x1+dx:dx:1000]; NX0 = length(X0);
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
% reflecteur horizontal
plot([x0 x1],[z0 z0],'r','LineWidth',2)
axis([0 1000 0 1000]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('Z (M)'), title('REFLECTEUR HORIZONTAL')


subplot(2,2,3), hold on, grid on
plot(0,0,'b','LineWidth',2), plot(0,0,'g',0,0,'b',0,0,'g')
legend('Réflexion','Diffraction','Temps réfléchis','Temps diffractés')
% diffractions
X = 0:10:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+z0^2)/V;
td1 = 2*sqrt((X-x1).^2+z0^2)/V;
plot(X,td0,'g',X,td1,'g')
X = 0:dx:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+z0^2)/V;
td1 = 2*sqrt((X-x1).^2+z0^2)/V;
for ix = 1:NX;
plot([X(ix) X(ix)],[0,td0(ix)],'g',[X(ix) X(ix)],[0,td1(ix)],'g')
end

% temps reflechis
X = x0:dx:x1; NX = length(X);
for ix = 1:NX
    plot(X(ix),0,'bo',X(ix),0,'bx',[X(ix) X(ix)],[0,t0],'b')
end
X0 = [0:dx:x0-dx , x1+dx:dx:1000]; NX0 = length(X0);
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
% reflexion horizontale
plot([x0 x1],[t0 t0],'b','LineWidth',2)
axis([0 1000 0 1]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('T (S)'), title('TEMPS REFLEXION - DIFFRACTION')

% pendage alpha
al = 30*pi/180; sal = sin(al); cal = cos(al); tal = tan(al); 

subplot(2,2,2), hold on, grid on
% rayons diffractes
X = 0:dx:1000; NX = length(X);
for ix = 1:NX
    plot([X(ix) x0],[0,x0*tal],'g',[X(ix) x1],[0,x1*tal],'g')
end
% rayons reflechis
X = x0+dx:dx:x1+2*dx; NX = length(X);
for ix = 1:NX
    plot(X(ix),0,'ro',X(ix),0,'rx',[X(ix) X(ix)*cal^2],[0,X(ix)*cal*sal],'r')
end
X0 = [0:dx:x0]; NX0 = length(X0);
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
% reflecteur pente
plot([x0 x1],[x0*tal x1*tal],'r','LineWidth',2)
axis([0 1000 0 1000]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('Z (M)'), title(['REFLECTEUR PENTE ' num2str(al/pi*180) ' °'])


subplot(2,2,4), hold on, grid on
% diffractions
zx0 = x0*tal; zx1 = x1*tal;
X = 0:10:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+zx0^2)/V;
td1 = 2*sqrt((X-x1).^2+zx1^2)/V;
plot(X,td0,'g',X,td1,'g')
X = 0:dx:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+zx0^2)/V;
td1 = 2*sqrt((X-x1).^2+zx1^2)/V;
for ix = 1:NX;
plot([X(ix) X(ix)],[0,td0(ix)],'g',[X(ix) X(ix)],[0,td1(ix)],'g')
end

% temps reflechis
X = x0+dx:dx:x1+2*dx; NX = length(X);
for ix = 1:NX
    plot(X(ix),0,'bo',X(ix),0,'bx',[X(ix) X(ix)],[0,2*X(ix)*sal/V],'b')
end
X0 = [0:dx:x0]; NX0 = length(X0);
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
% reflexion pentee
plot([x0/cal^2 x1/cal^2],[2*x0/cal^2*sal/V 2*x1/cal^2*sal/V],'b','LineWidth',2)
axis([0 1000 0 1]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('T (S)'), title('TEMPS REFLEXION - DIFFRACTION')


% reflecteur courbe
figure(2), clf 
% reflecteur concave
subplot(2,2,1), hold on, grid on
plot(0,0,'ko',0,0,'kx'), plot(0,0,'r','LineWidth',2), plot(0,0,'r',0,0,'g')
legend('Sources','Récepteurs','Réflecteur','Rayons réfléchis','Rayons diffractés',['V = ' num2str(V) ' M/S'])
% rayons diffractes
zx0 = sqrt(4*z0^2-(x0-500)^2)-z0;
zx1 = sqrt(4*z0^2-(x1-500)^2)-z0;
X = 0:dx:1000; NX = length(X);
for ix = 1:NX
    plot([X(ix) x0],[0,zx0],'g',[X(ix) x1],[0,zx1],'g')
end
% rayons reflechis
X = x0+dx:dx:x1-dx; NX = length(X);
for ix = 1:NX
    th = atan((X(ix)-500)/z0);
    plot(X(ix),0,'ro',X(ix),0,'rx',[X(ix) 2*z0*sin(th)+500],[0,z0*(2*cos(th)-1)],'r')
end
X0 = [0:dx:x0 x1:dx:1000]; NX0 = length(X0);
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
% reflecteur
X = x0:10:x1; NX = length(X);
plot(X,sqrt(4*z0^2-(X-500).^2)-z0,'r','LineWidth',2)
axis([0 1000 0 1000]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('Z (M)'), title('REFLECTEUR CONCAVE R = 1000 M')


subplot(2,2,3), hold on, grid on
plot(0,0,'b','LineWidth',2), plot(0,0,'g',0,0,'b',0,0,'g')
legend('Réflexion','Diffraction','Temps réfléchis','Temps diffractés')
% diffractions
X = 0:10:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+zx0^2)/V;
td1 = 2*sqrt((X-x1).^2+zx1^2)/V;
plot(X,td0,'g',X,td1,'g')
X = 0:dx:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+zx0^2)/V;
td1 = 2*sqrt((X-x1).^2+zx1^2)/V;
for ix = 1:NX;
plot([X(ix) X(ix)],[0,td0(ix)],'g',[X(ix) X(ix)],[0,td1(ix)],'g')
end
% temps reflechis
X = x0+dx:dx:x1-dx; NX = length(X);
for ix = 1:NX
    th = atan((X(ix)-500)/z0);
    plot(X(ix),0,'bo',X(ix),0,'bx',[X(ix) X(ix)],[0,2*z0*(2-1/cos(th))/V],'b')
end
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
% reflexion concave
X = x0:10:x1; NX = length(X);
th = asin((X-500)/z0/2);
plot(tan(th)*z0+500,2*z0*(2-1./cos(th))/V,'b','LineWidth',2)

axis([0 1000 0 1]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('T (S)'), title('TEMPS REFLEXION - DIFFRACTION')

% reflecteur convexe 
subplot(2,2,2), hold on, grid on
% rayons diffractes
zx0 = 3*z0-sqrt(4*z0^2-(x0-500)^2);
zx1 = 3*z0-sqrt(4*z0^2-(x1-500)^2);
X = 0:dx:1000; NX = length(X);
for ix = 1:NX
    plot([X(ix) x0],[0,zx0],'g',[X(ix) x1],[0,zx1],'g')
end
% rayons reflechis
X = x0-dx:dx:x1+dx; NX = length(X);
for ix = 1:NX
    th = atan((X(ix)-500)/z0/3);
    plot(X(ix),0,'ro',X(ix),0,'rx',[X(ix) 2*z0*sin(th)+500],[0,z0*(3-2*cos(th))],'r')
end
X0 = [0:dx:dx 1000-dx:dx:1000]; NX0 = length(X0);
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
axis([0 1000 0 1000]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('Z (M)'), title('REFLECTEUR CONVEXE R = 1000 M')
% reflecteur
X = x0:10:x1; NX = length(X);
plot(X,3*z0-sqrt(4*z0^2-(X-500).^2),'r','LineWidth',2)


subplot(2,2,4), hold on, grid on
% diffractions
X = 0:10:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+zx0^2)/V;
td1 = 2*sqrt((X-x1).^2+zx1^2)/V;
plot(X,td0,'g',X,td1,'g')
X = 0:dx:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+zx0^2)/V;
td1 = 2*sqrt((X-x1).^2+zx1^2)/V;
for ix = 1:NX;
plot([X(ix) X(ix)],[0,td0(ix)],'g',[X(ix) X(ix)],[0,td1(ix)],'g')
end
% temps reflechis
X = x0-dx:dx:x1+dx; NX = length(X);
for ix = 1:NX
    th = atan((X(ix)-500)/z0/3);
    plot(X(ix),0,'bo',X(ix),0,'bx',[X(ix) X(ix)],[0,2*z0*(3/cos(th)-2)/V],'b')
end
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
% reflexion convexe
X = x0:10:x1; NX = length(X);
th = asin((X-500)/z0/2);
plot(tan(th)*3*z0+500,2*z0*(3./cos(th)-2)/V,'b','LineWidth',2)
axis([0 1000 0 1]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('T (S)'), title('TEMPS REFLEXION - DIFFRACTION')

% courbure variable
figure(3), clf 
% reflecteur concave
subplot(2,3,1), hold on, grid on
% rayons diffractes
zx0 = sqrt(4*z0^2-(x0-500)^2)-z0;
zx1 = sqrt(4*z0^2-(x1-500)^2)-z0;
X = 0:dx:1000; NX = length(X);
for ix = 1:NX
    plot([X(ix) x0],[0,zx0],'g',[X(ix) x1],[0,zx1],'g')
end
% rayons reflechis
X = x0+dx:dx:x1-dx; NX = length(X);
for ix = 1:NX
    th = atan((X(ix)-500)/z0);
    plot(X(ix),0,'ro',X(ix),0,'rx',[X(ix) 2*z0*sin(th)+500],[0,z0*(2*cos(th)-1)],'r')
end
X0 = [0:dx:x0 x1:dx:1000]; NX0 = length(X0);
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
% reflecteur
X = x0:10:x1; NX = length(X);
plot(X,sqrt(4*z0^2-(X-500).^2)-z0,'r','LineWidth',2)
axis([0 1000 0 1000]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('Z (M)'), title(['RAYON COURBURE = ' num2str(2*z0) ' M'])


subplot(2,3,4), hold on, grid on
% diffractions
X = 0:10:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+zx0^2)/V;
td1 = 2*sqrt((X-x1).^2+zx1^2)/V;
plot(X,td0,'g',X,td1,'g')
X = 0:dx:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+zx0^2)/V;
td1 = 2*sqrt((X-x1).^2+zx1^2)/V;
for ix = 1:NX;
plot([X(ix) X(ix)],[0,td0(ix)],'g',[X(ix) X(ix)],[0,td1(ix)],'g')
end
% temps reflechis
X = x0+dx:dx:x1-dx; NX = length(X);
for ix = 1:NX
    th = atan((X(ix)-500)/z0);
    plot(X(ix),0,'bo',X(ix),0,'bx',[X(ix) X(ix)],[0,2*z0*(2-1/cos(th))/V],'b')
end
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
% reflexion concave
X = x0:10:x1; NX = length(X);
th = asin((X-500)/z0/2);
plot(tan(th)*z0+500,2*z0*(2-1./cos(th))/V,'b','LineWidth',2)

axis([0 1000 0 1]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('T (S)'), title('TEMPS REFLEXION - DIFFRACTION')


% centre en z=0 
subplot(2,3,2), hold on, grid on

% rayons diffractes
zx0 = sqrt(z0^2-(x0-500)^2);
zx1 = sqrt(z0^2-(x1-500)^2);
X = 0:dx:1000; NX = length(X);
for ix = 1:NX
    plot([X(ix) x0],[0,zx0],'g',[X(ix) x1],[0,zx1],'g')
end

% rayons reflechis
th = asin((500-x0)/z0);
plot(500,0,'ro',500,0,'rx',[500 500],[0 z0],'r',[500 x0],[0 z0*cos(th)],'r',[500 x1],[0 z0*cos(th)],'r')

X0 = [0:dx:500-dx 500+dx:dx:1000]; NX0 = length(X0);
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
% reflecteur
X = x0:10:x1; NX = length(X);
plot(X,sqrt(z0^2-(X-500).^2),'r','LineWidth',2)
axis([0 1000 0 1000]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('Z (M)'), title(['RAYON COURBURE = ' num2str(z0) ' M'])


subplot(2,3,5), hold on, grid on
% diffractions
X = 0:10:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+zx0^2)/V;
td1 = 2*sqrt((X-x1).^2+zx1^2)/V;
plot(X,td0,'g',X,td1,'g')
X = 0:dx:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+zx0^2)/V;
td1 = 2*sqrt((X-x1).^2+zx1^2)/V;
for ix = 1:NX;
plot([X(ix) X(ix)],[0,td0(ix)],'g',[X(ix) X(ix)],[0,td1(ix)],'g')
end
% temps reflechi
plot(500,0,'bo',500,0,'bx',[500 500],[0,t0],'b')
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
% reflexion
plot(500,t0,'b.','LineWidth',2)

axis([0 1000 0 1]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('T (S)'), title('TEMPS REFLEXION - DIFFRACTION')

% rayon courbure = z0/2
subplot(2,3,3), hold on, grid on
zc = z0/2; rc = z0-zc;
% rayons diffractes
zx0 = zc+sqrt(rc^2-(x0-500)^2);
zx1 = zc+sqrt(rc^2-(x1-500)^2);
X = 0:dx:1000; NX = length(X);
for ix = 1:NX
    plot([X(ix) x0],[0,zx0],'g',[X(ix) x1],[0,zx1],'g')
end
% rayons reflechis
X = x0-dx:dx:x1+dx; NX = length(X);
for ix = 1:NX
    th = atan((X(ix)-500)/zc);
    plot(X(ix),0,'ro',X(ix),0,'rx',[X(ix) -rc*sin(th)+500],[0 zc+rc*cos(th)],'r')
end
X0 = [0:dx:x0-2*dx x1+2*dx:dx:1000]; NX0 = length(X0);
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
% reflecteur
X = x0:10:x1; NX = length(X);
plot(X,sqrt(rc^2-(X-500).^2)+zc,'r','LineWidth',2)
axis([0 1000 0 1000]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('Z (M)'), title(['RAYON COURBURE = ' num2str(rc) ' M'])

subplot(2,3,6), hold on, grid on
% diffractions
X = 0:10:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+zx0^2)/V;
td1 = 2*sqrt((X-x1).^2+zx1^2)/V;
plot(X,td0,'g',X,td1,'g')
X = 0:dx:1000; NX = length(X);
td0 = 2*sqrt((X-x0).^2+zx0^2)/V;
td1 = 2*sqrt((X-x1).^2+zx1^2)/V;
for ix = 1:NX;
plot([X(ix) X(ix)],[0,td0(ix)],'g',[X(ix) X(ix)],[0,td1(ix)],'g')
end
% temps reflechis
X = x0-dx:dx:x1+dx; NX = length(X);
for ix = 1:NX
    th = atan((X(ix)-500)/zc);
    plot(X(ix),0,'bo',X(ix),0,'bx',[X(ix) X(ix)],[0,2*(zc./cos(th)+rc)/V],'b')
end
X0 = [0:dx:x0-2*dx x1+2*dx:dx:1000]; NX0 = length(X0);
for ix0 = 1:NX0
    plot(X0(ix0),0,'ko',X0(ix0),0,'kx')
end
% reflexion
X = x0:10:x1; NX = length(X);
th = asin((X-500)/rc);
plot(zc*tan(th)+500,2*(zc./cos(th)+rc)/V,'b','LineWidth',2)

axis([0 1000 0 1]), pbaspect('manual'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('T (S)'), title('TEMPS REFLEXION - DIFFRACTION')