介电结构的传递矩阵法

transfer matrix method for dielectric structures

我正在编写代码来模拟五层电介质结构的行为。使用的方法是TMM方法.. 输入是结构上的入射角。 对于这种结构,我期望的是临界角后的零透射率(此折射率堆栈为 65.1)和一段时间后出乎意料的峰值(这是由于光隧道效应)。 使用的公式取自论文:http://ieeexplore.ieee.org/document/6417944/

你知道为什么这个结构像法布里-珀罗结构吗?你认为我在循环中遗漏了什么吗?

非常感谢

`

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)*d(3)*sqrt((eps(4)-(n(1)^2)*(sin(theta)).^2));



m211=cos(delta1);
m212=j*sin(delta1)./nk2;
m221=j*nk2.*sin(delta1);
m222=cos(delta1);
%M2=[m211,m212;m221,m222];

m311=cos(delta2);
m312=j*sin(delta2)./nk3;
m321=j*nk3.*sin(delta2);
m322=cos(delta2);
%M3=[m311, m312;m321,m322];

m411=cos(delta3);
m412=j*sin(delta3)./nk4;
m421=j*nk4.*sin(delta3);
m422=cos(delta3);
%M4=[m411, m412;m421,m422];
Mtot=zeros(2,length(theta));

for i=1:length(theta)
M4=[m411(1,i),m412(1,i);m421(1,i),m422(1,i)];       
M3=[m311(1,i),m312(1,i);m321(1,i),m322(1,i)];
M2=[m211(1,i),m212(1,i);m221(1,i),m222(1,i)];
Mtot=M4*M3*M2;
T(i)=((nk5(i)/nk1(i))*(abs((2*nk1(i)/((Mtot(1,1)+Mtot(1,2).*nk5(i)).*nk1(i)+(Mtot(2,1)+Mtot(2,2).*nk5(i))))).^2));

end
figure(1);
plot(theta,T);
xlabel('incident angle');
ylabel('transmmission, s-pol');`

sin 以弧度为单位的角度作为参数,sind 以度为单位的角度。因此,您可以将所有 sin 更改为 sind(除了带有 delta 的那些)或将 theta 更改为弧度。后者最容易改变。像

theta=0:.01:pi/2;

就可以了,步子越小越好。同时将 plot 更改为

plot(theta*180/pi,T);

获得x-axis的学位。然后得到如下(for theta=0:.0001:pi/2;)