%OSEX05 calcule les hodochrones pour sujet controle continu ondes sismiques 2005
% voir : http://arnica.u-strasbg.fr/~jmm/Sismique/osex.html#C05
% fig1 = http://arnica.u-strasbg.fr/~jmm/Sismique/osex13.gif  
% fig2 = http://arnica.u-strasbg.fr/~jmm/Sismique/osex14.gif

clear all, clf

% vitesses, epaisseur, incidences critiques, gradient vitesse
V1 = 1500; H1 = 375;
V2 = 2000; H2 = 500; thc1 = asin(V1/V2);
V3 = 3000; thc2 = asin(V2/V3); thc21 = asin(V1/V3);
a = 2;

x = 0:1:4000;
x2 = x.^2;
% onde directe
td = x/V1;
% onde reflechie
tr = sqrt(4*H1^2+x2)/V1;
% onde conique interface 12
xc1 = 2*H1*tan(thc1), tc1 = 2*H1/cos(thc1)/V1,
x1 = x(find(x >= xc1));
t1 = x1/V2+2*H1*cos(thc1)/V1;
% onde conique interface 23
xc2 = 2*H2*tan(thc2)+2*H1*tan(thc21), tc2 = 2*H2/cos(thc2)/V2+2*H1/cos(thc21)/V1,
x2 = x(find(x >= xc2));
t2 = x2/V3+2*H2*cos(thc2)/V2+2*H1*cos(thc21)/V1;

figure(1), hold on
plot(x,td,'b',x,tr,'b',x1,t1,'r',x2,t2,'g','LineWidth',1.5), grid on
plot(xc1,tc1,'ro',xc2,tc2,'go')
xlabel('X EN M'), ylabel('T EN S'), title('HODOCHRONE')

% onde refractee dans couche a gradient vitesse lineaire
p = [1/V3:.00001:1/V2 , 1/V2];
th1 = asin(p*V1); th2 = asin(p*V2); R = 1/a./p;
x3 = 2*H1*tan(th1)+2*R.*cos(th2);
t3 = 2*H1./cos(th1)/V1-2/a*log(tan(th2/2));

figure(2), hold on
plot(x,td,'b',x,tr,'b',x3,t3,'r','LineWidth',1.5), grid on
plot(x3(1),t3(1),'ro',x3(length(p)),t3(length(p)),'ro')
xlabel('X EN M'), ylabel('T EN S'), title('HODOCHRONE')