%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Single Scattering Diffraction Stack Modeling
% Source on surface at (sx,0) and geophones at (gx,0)
% Homogeneous velocity model with velocity c and
% reflectivity model R(i,j)
%  c       - input- velocity of model in m/s
%  vel1    - input- velocity of heterogeneous model in m/s
% (nx,nz)  - input- Dimensions of Model (Make it square)
% (dx,dz)  - input- Spatial Sampling Intervals
% (dt,sx)  - input- Time sampling interval, source x-index 
% W(t)     -output- Ricker Wavelet with Central freq=fr
% D(x,t)   -output- Shot gather for a source at (sx,0)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
np=200;sx=1;sz=1;nx=100;nz=nx;
c=1000;dt=.005;dx=10;dz=dx;
fr=20;[rick]=ricker(np,fr,dt)
raymodel3;c=min(min(vel1));      % Compute traveltime tables traveltimesrc(:,:,isrc)
                                 % for velocity modle vel1;
nt=round(sqrt((dx*nx)^2+(dz*nz)^2)*2/c/dt+1)
nd2=round(nt)+1; W=zeros(2*nt,1);W(nt:nt+np-1)=rick(1:np);
x=[1:nx];z=[1:nz];

[X,Z]=meshgrid(x,z);
zz=(dz*Z).^2;X=(dx*X);                       % Compute X & Z Matrices
R=X*0;R(round(nz/2),round(nx/2))=1;          % Compute Refl. Model for pt scatterer
dvel1=(abs(diff(vel1,1)));mx=max(dvel1(:));  % Compute Refl Model for vel1;
dvel1=dvel1/mx;dvel1=(round(dvel1));imagesc(dvel1);colorbar
R=[dvel1;dvel1(nx-1,:)];                     % We lost a row in R by diff() so we add it back again

mig0a=R*0;icnt=0;m1b=mig0a;timer=round(traveltimesrc/dt)+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CREATE SHOT GATHER AND THEN MIGRATE IT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for isx=1:5:nx                                 % Loop Over Shots
icnt=icnt+1; TSX=timer(:,:,isx);               % Traveltime for a source at isx
  for  gx=1:nx;                                % Loop Over Traces
    gx
    xx=dx*gx-X;xx=xx.*xx;
    TXG=round(traveltimesrc(:,:,gx)/dt)+1;     % Traveltime for a src/geo at gx
    for t=1:nt;                                % Loop over time sample in a trace
      M=W(t-TSX-TXG+nd2).*R;
      D(gx,t)=sum(M(:));                       % Shot gather with source at isx
    end;
  end;

plot1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MIGRATE THE DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   gather1=diff(D(:,:),2,2);                   %Get CDP Shot Gather
   [m1b]=migvsp(timer,isx,nx,nz,dt,gather1);   %Migrate CDP Gather
   mig0a  = mig0a + m1b;
   imagesc(z*dx,x*dx,mig0a); xlabel('X (m)');ylabel('Y (m)');
   title(['Stacked Image for ',num2str(isx),' Shot Gathers']);pause(.25)
end
