function [D]=forw(R,nx,nz,nt,dx,dz,dt,c,sx,W,nd2,TGX)
% 
% Compute diffraction stack response D(gx,t) of 
% source at sx & reflectivity distribution given by R
%
x=[1:nx];z=[1:nz]; [X,Z]=meshgrid(x,z);
zz=(dz*Z).^2;X=(dx*X);             % Compute X & Z Matrices
D=zeros(nx,nt);cnst=1/c/dt;        %Define Constants
M=zeros(nx,nz);

for gx=1:nx;                       %Loop Over Traces
gx
for t=1:nt;                      %Loop over time sample/trace
M=W(t-TGX(:,:,gx)-TGX(:,:,sx)+nd2).*R;
D(gx,t)=sum(M(:));
  end;
end;
t=[0:nt-1]*dt;
imagesc(x*dx,t,D')
xlabel('X (m)');ylabel('Time (s)')

