%function image = rtm(xm,xs,xh,z,t,v,s1,wav,it0,npd);
function image = rtm1(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);
%     outape1(is,u);
% %     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);
%     outape2(is,u);
% %    ur(:,:,:,is)=u;
%     clear s2;
%     cpu2=cputime;
%     disp(['CPU time recv_extr:', num2str(cpu2-cpu1)]);
% end;

oxm=xm(1);nxm=length(xm);
dxm=xm(2)-xm(1);
oxs=xs(1);
dxs=xs(2)-xs(1);

image=zeros(nz,nxm,ns);

cpu1=cputime;
for is=1:ns;
%for is=31:40;
    load(['us_',num2str(is)]);
    load(['ur_',num2str(is)]);
 %   us1=fftshift(fft(us1,[],3));
 %   ur1=fftshift(fft(ur1,[],3));    
    img=zeros(nz,nh+2*npd);
    for it=1:nt;
        for ix=1:nh+2*npd;
            for iz=1:nz;
  %              img(iz,ix)=img(iz,ix)+conj(us1(iz,ix,it))*ur1(iz,ix,it);
                img(iz,ix)=img(iz,ix)+us1(iz,ix,nt-it+1)*ur1(iz,ix,nt-it+1);
            end;
        end;
    end;
    clear us1 ur1;

    ixm=floor(((is-1)*dxs+oxs+oxh-npd*dxh-oxm)/dxm)+1;
    
    image(:,ixm:ixm+nh-1+2*npd,is)=image(:,ixm-1:ixm+nh-2+2*npd,is)+img(:,:);
    disp(['shot #  ', num2str(is),' from ',num2str(ns),' xm(1) - xm(n):', num2str(ixm),' - ',num2str(ixm+nh-1+2*npd)]);
end;
cpu2=cputime;
disp(['CPU time fft(us and ur):', num2str(cpu2-cpu1)]);