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可知，每一個剖面的時間間隔就是
??????