%function image = rtm(xm,xs,xh,z,t,v,s1,wav,it0,npd);
function image = rtm(xm,xs,xh,z,t,v,wav,it0,npd);

% function img = rtm(xm,xs,xh,z,t,v,shot,wave);
%       xm --> x output domain
%       xs --> x shot positions
%       xh --> offset vector
%       v  --> must have (z,xm) dimension
%       ds and dxm must be equal

ns=length(xs);
nh=length(xh);
nm=length(xm);
nz=length(z);
nt=length(t);
[nzv,nxv]=size(v);

if nm~=nxv && nz~=nzv;
    error('Velocity and image domain have different sizes - Check for z and x');
end;

if ns~=1;
    dxs=xs(2)-xs(1);
else;
    dxs=0;
end;

oxh=xh(1);dxh=xh(2)-xh(1);

% for is=1:ns;
%     cpu1=cputime;
%     u = sht_down_2t8x(xs,is,xh,xm,z,t,v,wav,it0,npd);
%     us(:,:,:,is)=u;
%     cpu2=cputime;
%     disp(['CPU time shot_extr:', num2str(cpu2-cpu1)]);
% end;
 
for is=1:ns;
    cpu1=cputime;
    load(['shot_',num2str(is)]);
    s2=s1(1:nt,:);clear s1
    u = rec_down(xs,is,xh,xm,z,t,v,s2,npd);
    outape1(is,u);
%    ur(:,:,:,is)=u;
    clear s2;
    cpu2=cputime;
    disp(['CPU time recv_extr:', num2str(cpu2-cpu1)]);
end;

% cpu1=cputime;
% for is=1:is;
%     for ix=1:nh+2*npd;
%         for iz=1:nz;
%             us(iz,ix,:,is)=fftshift(fft(us(iz,ix,:,is)));
%             ur(iz,ix,:,is)=fftshift(fft(ur(iz,ix,:,is)));           
%         end;
%     end;
% end;
% cpu2=cputime;
% disp(['CPU time fft(us and ur):', num2str(cpu2-cpu1)]);

% % save us;
% % save ur;
% % 
% % load us;
% % load ur;

% cpu1=cputime;
% image = imaging(xs,xh,xm,us,ur,npd);
% %    size(img)
% %    ixm+nh-2+npd-(ixm-1-npd)
% %    (ixm+nh-1+npd)-(ixm-npd)
% cpu2=cputime;
% disp(['time imaging:', num2str(cpu2-cpu1)]);