matlab地震波数值模拟

%% 此程序为模拟各项同性弹性波波场（点源）
%% 时间上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）;
%*******************************************************************    
    
%***************************模型参数设置****************************%                
dx=5; dz=5;           %采样间隔
x0=0;z0=0;              %震源位置
t0=0.2;                %快照时间
l=length（t）;            %时间样点数目

x=-1000:dx:1000;
z=-1000:dz:1000;        %采样区间
m=length（x）;            
n=length（z）;            %采样点数
%% Setting Velocity & Density
xt=-1000:dx/2:1000+dx/2;
zt=-1000:dz/2:1000+dz/2;%采样区间
mt=length（xt）;  
nt=length（zt）;          %采样点数目

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）;
%**********************************************************************

%****************************稳定条件**********************************%
maxvp=max（max（vp））; maxvs=max（max（vs））;
if dt*sqrt（maxvp^2/dx^2+maxvs^2/dz^2）>1 
    error（‘参数设置有问题，差分格式不稳定‘）;
end
%**********************************************************************

%********************************数值计算******************************%
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;         %%时间差分参数
% 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;         %%时间上差分阶数
record1=zeros（l0-1n）;
record2=zeros（l0n）;

for k=1:l0;  %%空间上差分阶数
    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