%STOLT migre une coupe sismique a deport nul dans le domaine de Fourier
% voir : http://jmmeost.free.fr/Sismique/isis.html
% fig1 = http://jmmeost.free.fr/Sismique/stolt1.gif  

clear all, clf
colormap(pink)

% echantillonnage en t (s), f (Hz), x (m), cx (cycle/m), Nyquist
nt = 512; dt = .008;	tt = nt*dt;	 df = 1/tt; 	fn = 1/(2*dt); 	
nx = 64; dx = 50;	xx = nx*dx; 	dcx = 1/xx;	 cxn = 1/(2*dx);	
% vitesse et echantillonnage en z , cz, Nyquist
vit = 1000; 
nz = nt; dz = vit*dt; zz = vit*tt; 	dcz = 1/zz; czn = 1/(2*dz); 

% moitie superieure de la section
ttm = tt/2; itm = ttm/dt+1; 
zm = dz*(itm-1); czm = czn/2; izm = czm/dcz+1; fm = df*(izm-1); 

P = zeros(nt,nx); 

x = 0:dx:xx-dx; t = 0:dt:ttm ; [X,T] = meshgrid(x,t);

xx0 = xx/2; tt0 = tt/4; 
%fig1 onde plane lamdax=xx0 , periode=tt0 (enlever % sur ligne suivante)
P(1:itm,:) = cos(2*pi*(X/xx0-T/tt0));

%fig2 15 harmoniques, pas d'aliasing spatial (enlever % sur ligne suivante)
%nha = 15; for ic = 1:nha  P(1:itm,:) = P(1:itm,:) + cos(ic*2*pi*(X/xx0-T/tt0)); end, P = P/nha;

%fig3 30 harmoniques, aliasing spatial (enlever % sur ligne suivante)
%nha = 30; for ic = 1:nha  P(1:itm,:) = P(1:itm,:) + cos(ic*2*pi*(X/xx0-T/tt0)); end, P = P/nha;

%fig4 : impulsion en x=0 , t=.25s (enlever % sur ligne suivante)
%P(nt/16,1)=20;

%fig5 : impulsion en x=0 , t=1s (enlever % sur ligne suivante)
%P(nt/4,1)=20;

%fig6 : impulsion en x=1600m , t=1s (enlever % sur ligne suivante)
%P(nt/4,nx/2)=20; 

%dessine p(x,t)
subplot(2,2,1), pcolor(X,T,P(1:itm,:)) 
xlabel('X (M)'),ylabel('T (S)'), shading('flat')
set(gca,'TickDir','out', 'YDir', 'reverse') 
title('MIGRATION DE STOLT : P(X,T)')

%TF2D
P = fft2(P);

%dessine P(kx,omega)
subplot(2,2,2), axis([-cxn+dcx cxn 0 fm]) 
cx = -cxn+dcx:dcx:cxn; f = 0:df:fm; [CX,F] = meshgrid(cx,f); 
icx = [nx/2+2:nx 1:nx/2+1];
pcolor(CX,F,real(P(1:izm,icx)))
xlabel('KX (CYCLE/M)'), ylabel('FREQUENCE (HZ)'), shading('flat'), 
set(gca,'TickDir','out')
title('P(KX,OM) (partie réelle)')

%changement variable et interpolation lineaire
TFZ = 0.*P;
cz = dcz:dcz:czm; icx = 0;
%cas general -kxn < kx < kxn
for  cx = [0:dcx:cxn, -cxn+dcx:dcx:-dcx] 
   
%migration pour kx>kxn possible si toutes les pentes sont de meme signe
%voir : http://jmmeost.free.fr/Sismique/isex.html
%fig7 = http://jmmeost.free.fr/Sismique/stolt7.gif  
%for  cx =  [-2*cxn+dcx:dcx:-dcx]

   icx = icx+1;
   ff = vit*sqrt(cx^2 + cz.^2)/df+1; jf = fix(ff); q = (ff-jf).';   
   TFZ(2:izm,icx) = (1-q).*P(jf,icx) + q.*P(jf+1,icx);
   TFZ(nz:-1:nz-izm+2,icx) = (1-q).*P(nt-jf+2,icx) + q.*P(nt-jf+1,icx);
end

%dessine P(kx,kz) 
subplot(2,2,4), axis([-cxn+dcx cxn 0 czn])
cx = -cxn+dcx:dcx:cxn; cz = 0:dcz:czm; [CX,CZ] = meshgrid(cx,cz); 
icx = [nx/2+2:nx 1:nx/2+1];
pcolor(CX,CZ,real(TFZ(1:izm,icx)))
xlabel('KX (CYCLE/M)'), ylabel('KZ (CYCLE/M)'), shading('flat'), 
set(gca,'TickDir','out', 'DataAspectRatio',[1,1,1]) 
title('P(KX,KZ) (partie réelle)')

%TF2D inverse
P = real(ifft2(TFZ));

%dessine p(x,z)
subplot(2,2,3),  axis([0 xx 0 zm])
x = 0 :dx:xx-dx ; z = 0:dz:zm ; [X,Z] = meshgrid(x,z);
pcolor(X,Z,P(1:itm,:))
xlabel('X (M)'), ylabel('Z (M)'), shading('flat')
set(gca,'TickDir','out','YDir', 'reverse','DataAspectRatio',[1,1,1],'CLim',[0,1])
title('P(X,Z)   V = 1KM/S')