%?seismic?numerical?modeling
%?2D?time:?2?order?space:?8?order
%?no?boundary?condition?riker?wavelet
clc;
clear;
%%?many?parameters
%?grid?parameter?（m）?
dx=5.0;
dz=5.0;
nx=500;
nz=500;
x=（0:nx）*dx;
z=（0:nz）*dz;
%?time?parameter（s）
dt=0.001;??
nt=1000;
%?source?parameter?（m）
sx=1250;
sz=1250;
nsx=round（sx/dx）+1;
nsz=round（sz/dz）+1;
%seis_record=zeros（ntnx）;
%%?layers?and?velcoity（模型）
%??isotropic?medium
%v=3000*ones（nxnz）;?%單位是米每秒
v0=3000;
%?two?layers
%?v=3000*ones（nxnz）;
%?v（1:2001:end）=2500;
%?v（201:5001:end）=4000;
%?ricker?wavelet
f0=30;?%?（Hz）
t0=1.0/f0;
t=（1:nt）*dt;
Am=1000;??%大小待定
Wr=Am*（1-2*（pi*f0.*（t-t0））.^2）.*exp（-pi^2*f0^2.*（t-t0）.^2）;
p=zeros（nxnz）;
dpx=p;dpz=p;pold=p;pnew=p;
%?seis_record=zeros（ntnx）;
%?差分運算
for?ti=1:nt-2
????for?i=3:nx-3
????????for?j=3:nz-3
????????????dpx（ij）=（（4/3）*（p（i+1j）+p（i-1j））+（-1/12）*（p（i+2j）+p（i-2j））+...
????????????????（-5/2）*p（ij））/（dx^2）;
????????????dpz（ij）=（（4/3）*（p（ij+1）+p（ij-1））+（-1/12）*（p（ij+2）+p（ij-2））...
????????????????+（-5/2）*p（ij））/（d