%%?此程序為模擬各項同性彈性波波場（點源）
%%?時間上2階差分近似，空間上10階差分近似
clear;clc;
tic;

%*************************震源為Ricker子波**************************
dt=0.0001;t=0:dt:0.06;fp=30;
R=（1-2*（pi*fp*t）.^2）.*exp（-（pi*fp*t）.^2）;
%%plot（tR）;
%*******************************************************************????
????
%***************************模型參數(shù)設(shè)置****************************%????????????????
dx=5;?dz=5;???????????%采樣間隔
x0=0;z0=0;??????????????%震源位置
t0=0.2;????????????????%快照時間
l=length（t）;????????????%時間樣點數(shù)目

x=-1000:dx:1000;
z=-1000:dz:1000;????????%采樣區(qū)間
m=length（x）;????????????
n=length（z）;????????????%采樣點數(shù)
%%?Setting?Velocity?&?Density
xt=-1000:dx/2:1000+dx/2;
zt=-1000:dz/2:1000+dz/2;%采樣區(qū)間
mt=length（xt）;??
nt=length（zt）;??????????%采樣點數(shù)目

vp=zeros（mtnt）;????????%縱波速度
vs=zeros（mtnt）;????????%橫波速度
p=zeros（mtnt）;?????????%模型密度
%*****************************模型1************************************%
for?i=1:mt;
????for?j=1:nt;
????????vp（ij）=4300;
????????vs（ij）=2800;???%vp（ij）/sqrt（3）;
????????p（ij）=2.60;
????end
end
N?=?p.*vs.^2;?L?=?p.*（vp.^2-2*vs.^2）;
%**********************************************************************

%****************************穩(wěn)定條件**********************************%
maxvp=max（max（vp））;?maxvs=max（max（vs））;
if?dt*sqrt（maxvp^2/dx^2+maxvs^2/dz^2）>1?
????error（‘參數(shù)設(shè)置有問題，差分格式不穩(wěn)定‘）;
end
%**********************************************************************

%********************************數(shù)值計算******************************%
P=zeros（mn）;?Q=P;?S=P;?U=P;?W=P;
P1=P;?Q1=Q;?S1=S;?U1=U;?W1=W;
%?%C（1）=1.125;?C（2）=-0.04166667;?????????%%時間差分參數(shù)
%?C（1）=1.171875;?C（2）=-0.06510416;?C（3）=0.0046875;
C（1）=1.211243;?C（2）=-0.08972168;??C（3）=0.001384277;?C（4）=-0.00176566;?C（5）=0.0001186795;
m0=find（x==x0）;?n0=find（z==z0）;?l0=round（t0/dt）+1;
r=5;?????????%%時間上差分階數(shù)
record1=zeros（l0-1n）;
record2=zeros（l0n）;

for?k=1:l0;??%%空間上差分階數(shù)
????if?k<=l
????????P1（m0n0）=R（k）;?Q（m0n0）=R（k）;
????end
????for?i=r+1:m-r;
????????for?j=r+1:n-r;
????????????U（ij）=U1（ij）+C（1）*dt/p（2*i-12*j-1）*（1/dx*（P1（ij）-P1（i-1j））+1/dz*（S1（ij）-S1（ij-1）））+...
???????????????????????????C（2）*dt/p（2*i-12*j-1）*（1/dx*（P1（i+1j）-P1（i-2j））+1/dz*（S1（ij+1）-S1（ij-2）））+...
???????????????????????????C（3）*dt/p（2*i-12*j-1）*（1/dx*（P1（i+2j）-P1（i-3j））+1/dz*（S1（ij+2）-S1（ij-3）））+...
???????????????????????????C（4）*dt/p（2*i-12*j-1）*（1/dx*（P1（i+3j）-P1（i-4j））+1/dz*（S1（ij+3）-S1（ij-4）））+...
???????????????????????????C（5）*dt/p（2*i-12*j-1）*（1/dx*（P1（i+4j）-P1（i-5j））+1/dz*（S1（ij+4）-S1（ij-5）））;?
????????????W（ij）=W1（ij）+C（1）*dt/p（2*i2*j）*（1/dx*（S1（i+1j）-S1（ij））+1/dz*（Q1（ij+1）-Q1（ij）））+...
???????????????????????????C（2）*dt/p（2*i2*j）*（1/dx*（S1（i+2j）-S1（i-1j））+1/dz*（Q1（ij+2）+Q1（ij-1）））+...
???????????????????????????C（3）*dt/p（2*i2*j）*（1/dx*（S1（i+3