%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%?Single?Scattering?Diffraction?Stack?Modeling
%?Source?on?surface?at?（sx0）?and?geophones?at?（gx0）
%?Homogeneous?velocity?model?with?velocity?c?and
%?reflectivity?model?R（ij）
%??c???????-?input-?velocity?of?model?in?m/s
%??vel1????-?input-?velocity?of?heterogeneous?model?in?m/s
%?（nxnz）??-?input-?Dimensions?of?Model?（Make?it?square）
%?（dxdz）??-?input-?Spatial?Sampling?Intervals
%?（dtsx）??-?input-?Time?sampling?interval?source?x-index?
%?W（t）?????-output-?Ricker?Wavelet?with?Central?freq=fr
%?D（xt）???-output-?Shot?gather?for?a?source?at?（sx0）
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear?all;
np=200;sx=1;sz=1;nx=100;nz=nx;
c=1000;dt=.005;dx=10;dz=dx;
fr=20;[rick]=ricker（npfrdt）
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*nt1）;W（nt:nt+np-1）=rick（1:np）;
x=[1:nx];z=[1:nz];

[XZ]=meshgrid（xz）;
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（vel11）））;mx=max（dvel1（:））;??%?Compute?Refl?Model?for?vel1;
dvel1=dvel1/mx;dvel1=（roun