%DMO dessine l'operateur de Dip Move Out qui transforme les temps de 
% reflexion enregistres avec un demi-deport h le long de l'axe y des
% point-milieux en temps a deport nul quel que soit le pendage alpha
% voir : http://jmmeost.free.fr/Sismique/isis.html
% fig1 = http://jmmeost.free.fr/Sismique/nmo.gif
% fig2 = http://jmmeost.free.fr/Sismique/dmo.gif

clear all, clf

%echantillonage en point-milieux et deport
ym = 5000; y = 0:500:ym;  sy = size(y'); ny = sy(1);
hm = 2000; h = 100:100:hm; sh = size(h'); nh = sh(1);
dy0 = 100; tm = 4;
% vitesse et pendage alfa
vit = 1000; alfa = 0; %15*pi/180;

[Y,H] = meshgrid(y,h);

subplot(2,2,1), hold on, box on
title('TEMPS DE REFLEXION SUR PLAN PENDAGE ALPHA')
xlabel('Y(M)'), ylabel('H(M)'), zlabel('T(S)')
axis([0 ym 0 hm 0 tm]), set(gca, 'ZDir','reverse','CLim',[0 tm])

% calcul de t(y,h,alfa)
switch alfa
case 0
   T = 2/vit*sqrt(1000^2+H.^2); 
otherwise
   T = 2/vit*sqrt((Y+1000).^2*sin(alfa)^2+H.^2*cos(alfa)^2);
end
C = T; C(1:3,:) = nan*C(1:3,:);
% dessine t(y,h,alfa)
surf(Y,H,T,C), view([40,50])
switch alfa
case 0
   plot3(y,0*y,2/vit*1000*ones(sy)','r')
otherwise
   plot3(y,0*y,2/vit*(y+1000)*sin(alfa),'r') 
end

subplot(2,2,2), hold on, box on
title('Tn APRES CORRECTION DYNAMIQUE (NMO)')
xlabel('Y(M)'), ylabel('H(M)'), zlabel('T(S)') 
axis([0 ym 0 hm 0 tm]),  set(gca, 'ZDir','reverse', 'CLim',[0 tm])

% correction NMO
TN = real(sqrt(T.^2-(2/vit*H).^2)); 
C = TN; C(1:3,:) = nan*C(1:3,:);
% dessine tn(y,h,alfa)
surf(Y,H,TN,C), view([18,18])
switch alfa
case 0
   plot3(y,0*y,2/vit*1000*ones(sy)','r')
otherwise
   plot3(y,0*y,2/vit*(y+1000)*sin(alfa),'r')
end

subplot(2,2,4), hold on, box on
title('Tn POUR TOUS OFFSETS H (couleur=T(S))') 
xlabel('Y(M)'), ylabel('H(M)'), zlabel('T(S)')
axis([0 ym 0 hm 0 tm]) 
set(gca, 'ZDir','reverse', 'CLim',[0 tm]), colorbar,

% la correction NMO ne marche que pour alfa=0
surf(Y,H,TN,C), view([0,0])
switch alfa
case 0
   plot3(y,0*y,2/vit*1000*ones(sy)','r')
otherwise
   plot3(y,0*y,2/vit*(y+1000)*sin(alfa),'r')
end

subplot(2,2,3), hold on, box on
title('OPERATEURS DMO POUR 1 POINT-MILIEU Y'), 
xlabel('Y(M)'), ylabel('H(M)'), zlabel('T(S)')
axis([0 ym 0 hm 0 tm]), set(gca, 'ZDir','reverse', 'CLim',[0 tm]), 

% dessine t0(y,h=0,alfa)
view([174,16])
switch alfa
case 0
   surf(Y,H,2/vit*1000*ones(size(Y)))
otherwise
   surf(Y,H,2/vit*(Y+1000)*sin(alfa))
end

% operateur DMO
for iy = 7:7
   y0m = 2*hm^2/(vit*T(nh,iy)); 
   y0 = 0:dy0:y0m; y0 = [-max(y0):dy0:-dy0,y0]; [Y0,H0]=meshgrid(y0,h);
   switch alfa
   case 0
      T0 = 2/vit*1000 * real(sqrt(1-(Y0./H0).^2));
   otherwise
     T0 = 2/vit* real(sqrt(((y(iy)+1000).^2-H0.^2)*sin(alfa)^2)) .* real(sqrt(1-(Y0./H0).^2));
   end
   C0 = T0;
   for ih = 1:nh
      im = fix(2*h(ih)^2/(vit*T(ih,iy))/dy0); 
      i0 = fix(y0m/dy0); iyh = [1:i0-im-2,i0+im+3:2*i0+1];
      C0(ih,iyh) = nan*Y0(ih,iyh);
   end
   surf(Y0+y(iy),H0,T0,C0)
end