%入射速度波
a=0:pi/6:pi;
Vm=sin(a);
%線性漸變截面處各點高度
b(1)=0.4;b(9)=0.4;     
for i=1:1:7             
   b(i+1)=0.6;%等值變寬截面
end
%各單元之間的傳遞函數參數
Tu(1)=0.5;Td(1)=0.5;%存儲數據時多算兩組，以方便后面的遞推公式計算取值
Tu(10)=0.5;Td(10)=0.5;
for i=1:1:8
    Tu(i+1)=b(i+1)/(b(i+1)+b(i));
    Td(i+1)=b(i)/(b(i+1)+b(i));    
end
Tu(1)=0.5;Td(1)=0.5;%存儲數據時多算兩組，以方便后面的遞推公式計算取值
Tu(10)=0.5;Td(10)=0.5;
%第i單元第一時間初始條件
for i=1:1:31   %因為遞推公式中要用到樁底的信號,故這兒多一個初始值i=31
    Vu(i,1)=0; %j=1時,輸入信號為零,故各單元初始值為零
    Vd(i,1)=0;
end
%第i單元第j時間的速度
for j=2:1:70
    for i=2:1:30
       if i<=4
          H=[1 0]; h=[0 1];
        elseif i==5|i==6|i==12|i==13
          H=[2*Tu(i-3)*Tu(i-3)/Td(i-3) (Td(i-3)-Tu(i-3))*Td(i-4)/Tu(i-4)];
          h=[(Tu(i-3)-Td(i-3))*Tu(i-3)/Td(i-3) 2*Td(i-3)*Td(i-4)/Tu(i-4)];
       elseif  i==30      %樁底前一單元到樁底的傳遞函數
          H=[0 1];h=[0 0];
       else
          H=[1 0];h=[0 1];
       end
       Vu(i,j)=H*[Vu(i+1,j-1);Vd(i-1,j-1)]; 
       Vd(i,j)=h*[Vu(i+1,j-1);Vd(i-1,j-1)];
    end
   if j<=7
        Vd(1,j)=Vm(j); 
   elseif j>=10               %因為起始段為樁頂,時間j=0記為零時輸入信號時刻
        Vu(1,j)=Vu(2,j-1);%故從第一個變截面處反射到樁頂所需要的時間為j=10
        Vd(1,j+2)=Vu(1,j);
   else 
        Vd(1,j)=0;
   end
end
%樁頂接收到的速度波Vs
Vs=[];
for j=1:1:72
    if j<=7
      Vs(j)=Vd(1,j);
    else 
      Vs(j)=2*Vd(1,j);  
    end
end
j=1:1:72;
plot(j-1,Vs);%樁頂與第一個單元的時差為1個單位時間
grid on