声波方程数值模拟 FD方法

% 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