kirchhoff深度偏移

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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