有限差分正演模拟

nt=600;nx=301;nz=301;dx=10;dz=10;dt=0.001;
flag=6;
f=30;
tdlay=0.15;
sx=150;sz=150;
u0=zeros（nxnz）;u1=zeros（nxnz）;u2=zeros（nxnz）;v=zeros（nxnz）;


for i=1:nx
    for j=1:nz
        v（ij）=3000;
    end
end
%震源函数
s=zeros（nt1）;
%for
ts=0.3;
k=ts/dt;
it=（1:k）*dt;
    s=（1-2*（pi*f.*（it-1/f））.^2）.*exp（-（pi*f.*（it-1/f））.^2）;   
%end
pml=50;ix=0;iz=0;it=0;
px=0;pz=0;kk=0;izz=0;

nx2=nx+2*pml;nz2=nz+pml*2;
vel=zeros（nx2nz2）;
%定义速度模型
for ix=1:nx
    for iz=1:nz
        vel（ix+pmliz+pml）=v（ixiz）;
    end
end
for ix=1:pml
    for iz=1:nz2
        vel（ixiz）=vel（pml+1iz）;
        vel（ix+pml+nxiz）=vel（ix+pml+nx-1iz）;
    end
end
for ix=1:nx2
    for iz=1:pml
        vel（ixiz）=vel（ixpml+1）;
        vel（ixiz+nz+pml）=vel（ixiz+pml+nz-1）;
    end
end

P0=zeros（nx2nz2）;%前一时刻的波场
P1=zeros（nx2nz2）;%当前时刻的波场
P1t=zeros（nx2nz2）;%差分区域内部下一时刻的波场
P2t=zeros（nx2nz2）;%边界波场
wx=zeros（nx21）;
wz=zeros（nz21）;
w=zeros（nx2nz2）;
for ix=1:nx2
    wx（ix1）=1;
end
for iz=1:nz2
    wz（iz1）=1;
end
for ix=1:pml
    wx（ix1）=cos（pi/2*（pml-ix）/pml）;
    wx（nx2-ix1）=wx（ix1）;
end
for iz=1:pml
    wz（iz1）=cos（pi/2*（pml-iz）/pml）;
    wz（nz2-iz）=wz（iz1）;
end
for ix=1:nx2
    for iz=1:nz2
        w（ixiz）=wx（ix）*wz（iz）;
    end
end


isnap=30;


a=zeros（81）;
c=zeros（88）;
a（11）=-2;%2jie
a（21）=-2.5;%4jie
a（31）=-2.72222;%6jie
a（41）=-2.8472;%8jie
a（51）=-2.9272;%10jie
a（61）=-5369/1800;%12jie
a（71）=-266681/88200;%14jie
a（81）=-1077749/352800;%16jie
c（11）=1;%2jie
c（21）=4/3;c（22）=-1/12;%4jie
c（31）=3/2;c（32）=-3/20;c（33）=1/90;%6jie
c（41）=8/5;c（42）=-1/5;c（43）=8/315;c（44）=-1/560;%8jie
c（51）=5/3;c（52）=-5/21;c（53）=5/126;c（54）=-5/1008;c（55）=1/3150;%10jie
c（61）=12/7;c（62）=-15/56;c（63）=10/189;c（64）=-1/122;c（65）=2/1925;c（66）=-1/16632;%12jie
c（71）=7/4;c（72）=-7/24;c（73）=7/108;c（74）=-7/528;c（75）=7/3300;c（76）=-7/30888;c（77）=1/84084;%14jie
c（