function image=imaging(xs,xh,xm,us,ur,npd);

% function img=imaging(us,ur,nh,npd)
%
% us,ur --> nz,nx,nf,ns
% 
nz=size(us,1);
nh=size(us,2);
nw=size(us,3);
ns=size(us,4);

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

image=zeros(nz,nxm);
for is=1:ns;
    img=zeros(nz,nh);
    for iw=1:nw;
        for ix=1:nh;
            for iz=1:nz;
                img(iz,ix)=img(iz,ix)+conj(us(iz,ix,iw,is))*ur(iz,ix,iw,is);
            end;
        end;
    end;

    ixm=floor(((is-1)*dxs+oxs+oxh-npd*dxh-oxm)/dxm)+1;
    
    image(:,ixm:ixm+nh-1)=image(:,ixm-1:ixm+nh-2)+img(:,:);
    disp(['shot #  ', num2str(is),' from ',num2str(ns)]);
end;
image=real(image);