用传输矩阵方法编写的matalb程序

na=2.3;nb=1;nc=1.38;n1=1;n2=1;%边界为空气
for d=400:1600;   %波长范围
c3=0;  % 第一层入射角
c1=asin(n1*sin(c3)/na);%第一层折射角，第二层入射角
c2=asin(na*sin(c1)/nb);%第二层折射角
%c4=asin(nb*sin(c2)/n2);%出射角
c5=asin(nb*sin(c2)/nc);
d1=1064;%中心波长
a=d1/(4*na);b=d1/(4*nb);c=2*d1/(4*nc);%每层介质的光学厚度均为某一参考波长的1/4,a,b,c介质层的几何厚度
Ba=2*pi*na*a*cos(c1)/d;%特征拒阵里的项
Bb=2*pi*nb*b*cos(c2)/d;%特征拒阵里的项
Bc=2*pi*nc*c*cos(c5)/d;
f=4*pi*1e-7;       %真空磁导率
e=1e-9/(36*pi);    % 真空介电常数
m=sqrt(e/f);        %真空光速
za=m*cos(c1)*na;zb=m*cos(c2)*nb;
z1=f*cos(c3)*n1;z2=f*cos(c5)*n2;
p1=cos(Bb);p2=-i*sin(Bb)/zb;
p3=-i*zb*sin(Bb);p4=cos(Bb);
P=[p1 p2;p3 p4];
y1=cos(Bc);y2=-i*sin(Bc)/(m*nc*cos(c5));y3=-i*m*nc*cos(c5)*sin(Bc);y4=cos(Bc);
Y=[y1 y2;y3 y4];
q1=cos(Ba);
q2=-i*sin(Ba)/za;
q3=-i*za*sin(Ba);
q4=cos(Ba);
Q=[q1 q2;q3 q4];
O=Q*P;%OC=P*Q;
O1=O^8*Y*O^8;
O11=O1(1,1);O12=O1(1,2);O13=O1(2,1);O14=O1(2,2);
z1=sqrt(e/f)*n1*cos(c3);z2=sqrt(e/f)*n2*cos(c5);
t=2*z1/(z1*(O11+z2*O12)+O13+z2*O14);%透射
r=(z1*(O11+z2*O12)-O13-z2*O14)/(z1*(O11+z2*O12)+O13+z2*O14);%反射
t1=abs(t);%绝对值
r1=abs(r);
s=d-399;
k(1,s)=t1
x(1,s)=r1
end
d=400:1600;
%subplot(1,2,1)
plot(d,k);xlabel('Incident wavelength(nm)','fontsize',19);
ylabel('Transmission','fontsize',19);
gtext('na=2.3,nb=1,nc=1.38,N=8,m=2','fontsize',14)
%subplot(1,2,2)
%plot(d,x)