%VRMS dessine la vitesse RMS correspondant a V(z)=V0+az et l'erreur sur
% le temps de parcours entre une source en ZS et un recepteur en z = 0
% due a l'approximation hyperbolique avec VRMS
% voir : http://jmmeost.free.fr/Sismique/osis.html
% fig = http://jmmeost.free.fr/Sismique/vrms.gif

clear all; clf;
c = ['brgmc']; co =['ob';'or';'og';'om';'oc'];

V0 = 1000; a = 1;
XM = 5000; ZM = 5000;

z = 1:1:ZM; V = V0+a*z;

subplot(1,2,1), hold on, grid on
title('VRMS et VM POUR V(Z) = 1000+Z M/S')
xlabel('VITESSE (M/S)'), ylabel('PROFONDEUR (M)')
set(gca,'YDir','reverse','DataAspectRatio',[1,1,1])

% VRMS et vitesse moyenne
VR = sqrt(a*(V0*z+a*z.^2/2)./log(1+a*z/V0)); 
VM = a*z./log(1+a*z/V0);

plot(V,z,'b',VR,z,'r',VM,z,'g')

subplot(1,2,2), hold on, grid on
title('ERREUR SUR TEMPS DE PARCOURS DE ZS A Z = 0')
ylabel('ERREUR RELATIVE MAX (%) ,  T(VRMS) - T(V0+AZ) (MS)'), xlabel('DISTANCE HORIZONTALE (M)')
axis([0,XM,0,100])

% source en ZS
for ZS = 1000:1000:5000
VS = V0+a*ZS;
VRS2 = a*(V0*ZS+a*ZS^2/2)/log(1+a*ZS/V0);

subplot(1,2,1), plot(sqrt(VRS2),ZS,'ro');

% rayons montants d'incidence ths en ZS, th0 en z = 0
ths = [0.1:.1:pi/2,pi/2];
p = sin(ths)/VS; R = 1./p/a;
th0 = asin(p*V0); X0 = R.*(cos(th0)-cos(ths)); 
T0 = log(tan(ths/2)./tan(th0/2))/a*1000;

% approximation hyperbolique avec VRMS
t0 = log(VS/V0)/a ; t2 = t0^2; X2 = X0.*X0; 
TR = sqrt(t2+X2/VRS2)*1000;

% dessine l'ecart de temps 
subplot(1,2,2)
plot([0,X0],[0,TR-T0],c(ZS/1000))
end

%erreur relative maximum
for ZS = 1000:1000:5000
VS = V0+a*ZS;
VRS2 = a*(V0*ZS+a*ZS^2/2)/log(1+a*ZS/V0);
th0 = asin(V0/VS); X0 = VS*(cos(th0))/a;
T0 = -log(tan(th0/2))/a*1000;
t0 = log(VS/V0)/a ; t2 = t0^2; X2 = X0*X0;
TR = sqrt(t2+X2/VRS2)*1000;
plot(X0,100*(TR-T0)/T0,co(ZS/1000,:))
end

subplot(1,2,1), legend('V0+AZ','VRMS','VM')
subplot(1,2,2), legend('ZS = 1 KM','2 KM','3 KM','4 KM','5 KM','% ERREUR')