transfer matrix method for dielectric structures

Question

I am writing a code to simulate the behaviour of a five-layer dielectric structure. The method used is the TMM methods.. Inputs are the incident angle on the structure. With this strcture, what I expected is a zero transmission after the critical angle (which is 65.1 for this refractive index stacks) and an unexpected peak after a while (this is due to an optical tunneling effect). The formula used are taken from paper: http://ieeexplore.ieee.org/document/6417944/

Do you know why the structure acts like a Fabry-Perot instead? Am I missing something in the loop in your opinion?

Thank you very much

`

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');`

Answer 1

sin takes as argument an angle in radians, sind takes the angle in degrees. So you can either change all sin to sind (expect the ones with delta) or change theta to radian. The latter is easiest to change. Something like

theta=0:.01:pi/2;

will be good enough, smaller steps are of course better. Also change the plot to

plot(theta*180/pi,T);

to have the x-axis in degree. Then the following is obtained (for theta=0:.0001:pi/2;)