四层声波波场有限差分程序

function yxcf
% …………………………此程序用于模拟声波波场…………………………
clear;clc;
% 雷克子波
figure（1）
dt=0.0001;t=（0:150）*dt;f=30; %这里从0开始就是半个雷克子波，t是41个值，f表示主频为30Hz。这里该波的时间长度是0.04秒
R=（1-2*（pi*f*t）.^2）.*exp（-（pi*f*t）.^2）; %R也是41个值，R乘以一个数可以改变最大振幅的大小
% plot（tR）
% title（‘子波模型‘）

% …………………………模型参数设置…………………………
dx=2;dy=2;%采样间隔
x=-100:dx:100;y=-100:dy:100;%采样区间
x0=0;y0=0;t0=0.08;%炮点激发地点和激发时刻
m=length（x）; %x数组的长度
n=length（y）; %y数组的长度，确定整个计算区域为101点见方的区域。

%…………………………数值计算…………………………
data0=zeros（mn）; %data0表示i-1时刻的波场值，用m和n表示矩阵行列
data1=data0; %data1表示i时刻的波场值
data2=data0; %data2表示i+1时刻的波场值
v=zeros（mn）;
m0=find（x==x0）;n0=find（y==y0）;h0=round（t0/dt）+1; %（m0n0）是坐标原点（00）在x和y矩阵中位置的标号m0=n0=51
                                                 %（利用在数组中分别找出x=0和y=0的下标）。h0是801
%…………………………水平四层模型…………………………
s1=find（y==-50）;
s2=find（y==0）;
s3=find（y==50）;
for i=1:m
    for j=1:n
        if i<=s1
            v（ij）=700;
        else if i<=s2
                v（ij）=1000;
            else if i<=s3
                    v（ij）=1400;
                else if i>s3
                        v（ij）=2800;
                    end
                end
            end
        end
    end
end
%………………………………………………………………
global T
T=2000;%设定采样时间单位ms
for k=1:T      %k=1:h0这里k表示的是第几个时间剖面，k=1时表示起震，
                 %data1（m0n0）=R（k）这句话表示把雷克子波的最大值给（0，0）点，这是起震点，后面随着
                 %k的值增大，波就一点一点慢慢向前推进。且由R（k）-R（k-1）和dt=0.001可知，每一个剖面的时间间隔就是