基于TMM的五层介质结构模拟异常现象技术问询

我正在用**传输矩阵法（TMM）**编写代码，模拟一个五层介质结构的光学特性，输入变量是入射角。根据这个折射率堆叠的参数，理论临界角是65.1度——我预期临界角之后透射率会降到0，后续会出现光学隧穿效应导致的异常峰。公式参考了一篇IEEE学术论文，但目前的模拟结果却显示这个结构表现得像法布里-珀罗腔，完全不符合预期。

我怀疑是不是代码的循环部分遗漏了什么关键步骤？麻烦各位帮忙看看我的代码，非常感谢！

代码片段： #

N=5; %numbers of layer
theta=0:90;
lambda=800;%wavelength in vacuum in nm
eps1=2.28;eps2=1.87;eps3=2.28;eps4=1.87;eps5=2.28;
d=[600,400,600];%thickness of layer in nm for layer 2,3,4
mu=[1,1,1,1,1]; %permeability of every layer
eps=[eps1,eps2,eps3,eps4,eps5];
n=zeros(1,5);%refractive index
T=zeros(1,length(theta));
for s=1:5
    n(s)=sqrt(eps(s)*mu(s));
end
nk1=sqrt((eps(1)-(n(1)^2)*(sin(theta)).^2));%optical admittance inc layer
nk2=sqrt((eps(2)-(n(1)^2)*(sin(theta)).^2));%optical admittance second layer
nk3=sqrt((eps(3)-(n(1)^2)*(sin(theta)).^2));%optical admittance third layer
nk4=sqrt((eps(4)-(n(1)^2)*(sin(theta)).^2));%optical admittance fourth layer
nk5=sqrt((eps(5)-(n(1)^2)*(sin(theta)).^2));%optical admittance fith layer
delta1=((2*pi)/lambda)*d(1)*sqrt((eps(2)-(n(1)^2)*(sin(theta)).^2));
delta2=((2*pi)/lambda)*d(2)*sqrt((eps(3)-(n(1)^2)*(sin(theta)).^2));
delta3=((2*pi)/lambda)...

内容的提问来源于stack exchange，提问作者user3050442