把兩個源代碼放在一起，運行即可了！
這是
波場計算的源代碼
function seismography_disp=wy_vsp(p)
% set the wavelet_length as a odd number
%wavelet_time=0.1;
wavelet=p.wavelet;
wavelet=wavelet.*2;
fs=p.fs;
density=p.density;
Q=p.Q;
v=p.v;
time=p.time;
dh=p.dh;
layer_number=p.layer_number;
R_n=layer_number-1;
wavelet_length=p.wavelet_length;
%計算 區間復數速
for j=1:layer_number
    vc(j,:)=v(j,:)./(1-i/(2*Q(j)));
end
%計算反射系數
R=zeros(layer_number,wavelet_length);
for j=1:R_n
     R(j,:)=(density(j+1).*real(vc(j+1,:))-density(j).*real(vc(j,:)))./(density(j+1).*real(vc(j+1,:))+density(j).*real(vc(j,:)));
%      if max(R(j,:)<=0.00001
%          R(j,:)=0;
%      end
end
% 計算傳播矩陣
A1=zeros(R_n,wavelet_length);
A2=zeros(R_n,wavelet_length);
A3=zeros(R_n,wavelet_length);
A4=zeros(R_n,wavelet_length);
temp = 1:wavelet_length;
w=double(temp-1)*2*pi/time;
q=zeros(layer_number,wavelet_length);
for j=1:layer_number
    q(j,:)=exp( -w.*dh(j)./(2*Q(j).*v(j,temp)));%Q
end
for j=2:R_n
    t2=exp( -i*w.*dh(j)./v(j,temp));%up
    t1=exp( -i*w.*dh(j)./v(j,temp));% down
    r1=1./(1+abs(R(j,temp)));%defraction
    r2=R(j,temp)./(1+abs(R(j,temp)));%reflection
    A1(j,temp)=r1.*q(j+1).*t1;
    A2(j,temp)=r2.*q(j+1).*t2;
    A3(j,temp)=r2.*q(j).*t1;
    A4(j,temp)=r1.*q(j).*t2;
end
clear R 
%define ups and downs
up=zeros(R_n,wavelet_length);
down=zeros(R_n,wavelet_length);
tempup=zeros(R_n,wavelet_length);
tempdown=zeros(R_n,wavelet_length);
seismography=zeros(R_n,wavelet_length);
downgraphy=zeros(R_n,wavelet_length);
upgraphy=zeros(R_n,wavelet_length);
down(2,:)=wavelet;
seismography(2,:)=wavelet;
for j=1:80  % 計算傳播的次數
    %up(1,:)=0;down(R_n,:)=0;%
    t=1;
    tempdown(t+1,:)=exp( -w.*dh(1)./(2*Q(1).*v(1,temp))) .* exp( -i*w.*dh(1)./v(1,temp)).* up(t,:);
    for t=2:R_n-1
        tempup(t-1,:)=A1(t,:).*up(t,:)+A2(t,:).*down(t,:);
        tempdown(t+1,:)=A3(t,:).*up(t,:)+A4(t,:).*down(t,:);
    end
    t=R_n;
    tempup(t-1,:)=A1(t,:).*up(t,:)+A2(t,:).*down(t,:);
    up=tempup;
    down=tempdown;
%     figure;imagesc (abs(tempup)); figure(gcf);title('up');
%     figure;imagesc (abs(tempdown)); figure(gcf);title('down');
    seismography=seismography+up+down;
    upgraphy=upgraphy+up;
    downgraphy=downgraphy+down;
end
clear tempdown tempup down up 
% transform the wavelets to time domain
wavel=fs*time-1;
seismography_disp=zeros(R_n,wavel);
upgraphy_disp=zeros(R_n,wavel);
downgraphy_disp=zeros(R_n,wavel);
for j=1:R_n
    seismography_disp(j,1)=seismography(j,1);
    upgraphy_disp(j,1)=upgraphy(j,1);
    downgraphy_disp(j,1)=downgraphy(j,1);
    for jj=2:wavelet_length
        seismography_disp(j,jj)=seismography(j,jj);
        upgraphy_disp(j,jj)=upgraphy(j,jj);
        downgraphy_disp(j,jj)=downgraphy(j,jj);
        seismography_disp(j,wavel-jj+2)=seismography_disp(j,jj)';
        upgraphy_disp(j,wavel-jj+2)=upgraphy_disp(j,jj)';
        downgraphy_disp(j,wavel-jj+2)=downgraphy_disp(j,jj)';
    end
    seismography_disp(j,:)=ifft(seismography_disp(j,:));
    upgraphy_disp(j,:)=ifft(upgraphy_disp(j,:));
    downgraphy_disp(j,:)=ifft(downgraphy_disp(j,:));
end
showtrace(seismography_disp);
showtrace(upgraphy_disp);
showtrace(downgraphy_disp);
 
這是 參數設置的源代碼：
function p=para()
time=5;
p.time=time;
[p.wavelet,p.fs]=leikewavelet(50);
fs=p.fs;
wavelet=p.wavelet;
% wavelet
wavelet_length=length(p.wavelet);
if int32(wavelet_length/2)*2==wavelet_length
    p.wavelet(wavelet_length+1)=0;
end
wavelet_initial=zeros(1,time*fs);
wavelet_initial(1:wavelet_length)=wavelet;
fft_wavelet=fft(wavelet_initial);
[temp,wavelet_length]=max(abs(fft_wavelet));
clear temp
wavelet_length=wavelet_length*3;
wavelet_length=int32(wavelet_length);
%fft the wavelet and use the native part of the fft_wavelet as the initial
%wavelet in the frequency domain.
wavelet_initial = fft_wavelet(1:wavelet_length); % get the initial wavelet
wavelet_length = length(wavelet_initial);
p.wavelet_length = wavelet_length;
p.layer_number=50;
p.v=zeros(p.layer_number,p.wavelet_length);
p.density=zeros(p.layer_number,1);
% layer's interval
dh=zeros(1,p.layer_number);
p.dh=dh+100;
p.wavelet=wavelet_initial;
% %計算各個層速度，密度,傳播時間time
% vv=2200;% constant velocity;
% pp=2000;%
% for j=1:p.layer_number
%     p.v(j,:)=vv+(j-1)*10;%linear additve velocity
%     p.density(j)=pp+(j-1)*2;%
% end
%p.Q=200;% Quanity ceofficient
%%%%%設計了5個地層,velocity and density
p.v(1:4,:)=1500;density(1:4)=2000;p.Q(1:4)=200;
p.v(5:20,:)=3000;density(5:20)=3000;p.Q(5:20)=500;
p.v(21:25,:)=2000;density(21:25)=1500;p.Q(21:25)=100;
p.v(26:35,:)=3000;density(26:35)=3000;p.Q(26:35)=300;
p.v(36:50,:)=5000;density(36:50)=5000;p.Q(36:50)=500;
p.density=density;
%%%%%%% 速度隨圓頻率變化
w=1:p.wavelet_length;
for j=1:p.layer_number
    p.v(j,w)=p.v(j,w).*(1+0.00001.*w);
end