ViscoAcoustic_2D_FDTD_SG2.m

clc,clear
% ViscoAcoustic wave simulation using finite difference. with staggered-grid.
% By zhaoqingwei
% Chengdu University of Technology (CDUT), 2021-2025
%%%%%%%%%%%%%%%%
% P--Vx--P
% |      |
% Vz-----vz
% |      |
% P--Vx--P
%%%%%%%%%%%%%%%

FM=20;
DT=0.001;dt=DT;
T=2;
nx=600;
nz=600;
DH=10;dx=DH;dz=DH;

s=wavelet(FM,DT,T)*1e5;
vx0=zeros(nz,nx);vx1=vx0;
vz0=vx0;vz1=vx0;
p0=vx0;p1=vx0;
v=ones(nz,nx)*4000;Q=ones(nz,nx)*100;vq=v./sqrt(Q*FM);
z0=round(nz/2);x0=round(nx/2);

nn=3;
dxd=FDcoeffDx(nn);ddz0=dxd';ddz1=[dxd 0]';ddx0=ddz0';ddx1=ddz1';
TX=dt/dx;TZ=dt/dz;
for	t=DT:DT:T
    disp(t);
    k=round(t/DT);
    vz1=vz0+TZ*imfilter(p0,ddz0);
    vx1=vx0+TX*imfilter(p0,ddx0);
    vzt=TZ*imfilter(p0,ddz0)/dt;
    vxt=TX*imfilter(p0,ddx0)/dt;
    
    p1=p0+v.*v.*(TZ.*imfilter(vz1,ddz1)+TX.*imfilter(vx1,ddx1))+vq.*vq.*(TZ.*imfilter(vzt,ddz1)+TX.*imfilter(vxt,ddx1));
    
    p1(z0,x0)=p1(z0,x0)+s(k);
    p0=p1;
    vz0=vz1;vx0=vx1;
    
    if mod(t,0.5)==0
        figure();
        imagesc(real(p1));
    end
end