clear;
T=4 ; %輸入介質的層數
n=[1,2.17,2.89,3.42] ;%輸入每層介質的折射率
Angle1=0; % 入射角
Depth=[0,0.15,7,3];%輸入每層介質的長度
WaveNumber=500:1:5000;
j=1; %坐標刻度間隔
for t=1:T                     
Angle(t)=asin(n(1)*sin(Angle1)/n(t));    
end
l=1;
Ds0=[1,1;n(1)*cos(Angle(1)),-n(1)*cos(Angle(1))];
Dp0=[cos(Angle(1)),cos(Angle(1));n(1),-n(1)];
DsT=[1,1;n(T)*cos(Angle(T)),-n(T)*cos(Angle(T))];
DpT=[cos(Angle(T)),cos(Angle(T));n(T),-n(T)];
 for wavenumber=WaveNumber(1):j:WaveNumber(length(WaveNumber))
    wavelength=10000/wavenumber;
    DS=1;
    DP=1;
    for k=2:T-1              
    Ds=[1,1;n(k)*cos(Angle(k)),-n(k)*cos(Angle(k))];
    Dp=[cos(Angle(k)),cos(Angle(k));n(k),-n(k)];
    Phasic(k)=2*pi*n(k)*Depth(k)*cos(Angle(k))/(wavelength);
    Pl=[exp(i*Phasic(k)),0;0,exp(-i*Phasic(k))];
    DSn=Ds*Pl*inv(Ds);
    DS=DS*DSn;
    DPn=Dp*Pl*inv(Dp);
    DP=DP*DPn;
    end      
   DSFINAL=inv(Ds0)*DS*DsT;
   DPFINAL=inv(Dp0)*DP*DpT;
   MS11=DSFINAL(1,1);
   MS21=DSFINAL(2,1);
   MP11=DPFINAL(1,1);
   MP21=DPFINAL(2,1);
   rs=MS21/MS11;
   rp=MP21/MP11;
   RS=abs(rs)^2;
   RP=abs(rp)^2;
   R(l)=(RS+RP)/2;
   l=l+1; 
end
 plot(WaveNumber,R);