clear
close all
maxw=2.5e15;
minw=1e11;
Nw=2000;
w=linspace(minw,maxw,Nw);
kj=1;
for j=1:Nw;

lam0=2*pi*3e14./w(j); % vacuum wavelength in micron

%% output and input materials
nout=1.5;
nin=1;
k0out=2*pi*nout./lam0;
k0in=2*pi*nin./lam0;

Np=20;

thetain=pi/4;
ky=k0in*sin(thetain);
thetaout=asin(ky./k0out)
%% unit cell  (AB)%%%%%%%%

nA=2;             % refractive index of material A
dA=0.25;            % thickness of layer A in micron
k0A=2*pi*nA./lam0;  % modulus of the wavevector in layer A ( micron^(-1))

nB=1.428;     % refractive index of material B
dB=0.35;    % thickness of layer B in micron
k0B=2*pi*nB./lam0;  % modulus of the wavevector in layer B ( micron^(-1))

a=dA+dB;    % unit cell thickness




kxA=sqrt(k0A^2-ky^2);   % x component of the wavevector in layer A ( micron^(-1))
kxB=sqrt(k0B^2-ky^2);   % x component of the wavevector in layer B (micron^(-1))
kxin=sqrt(k0in^2-ky^2);
kxout=sqrt(k0out^2-ky^2);

%% Transfer matrix of the unit cell  %%%%
M=zeros(2,2);

PB= [ exp(-i.*kxB.*dB)   0;
      0                  exp(i.*kxB.*dB)];
 
IBA=0.5*[ 1+(kxB/kxA)   1-(kxB/kxA);
      1-(kxB/kxA)   1+(kxB/kxA)];
  
PA= [ exp(-i.*kxA.*dA)   0;
      0                  exp(i.*kxA.*dA)];
 
IAB=0.5*[ 1+(kxA/kxB)   1-(kxA/kxB);
      1-(kxA/kxB)   1+(kxA/kxB)];
  
M=IAB*PA*IBA*PB;

%% finite structure
Iout=0.5*[ 1+(kxout/kxB)   1-(kxout/kxB);
      1-(kxout/kxB)   1+(kxout/kxB)];

  Iout1=0.5*[ 1+(kxout/kxA)   1-(kxout/kxA);
      1-(kxout/kxA)   1+(kxout/kxA)];
Iin=0.5*[ 1+(kxA/kxin)   1-(kxA/kxin);
      1-(kxA/kxin)   1+(kxA/kxin)];
  Iin1=0.5*[ 1+(kxB/kxin)   1-(kxB/kxin);
      1-(kxB/kxin)   1+(kxB/kxin)];
  
  NM=(Iin*PA*IBA*PB*(M^(Np-1))*Iout);
  %NM=M^Np;
%% eigen values %%
Ka=acos(0.5*(M(1,1)+M(2,2)))./a;
%K=0;
%KY=0;
%if abs(imag(Ka))<2;
tt(kj)=0;
    if abs(real(Ka))>0;
        if abs(real(Ka))<pi./a;
K(kj)=real(Ka);
freq(kj)=w(j);
KY(kj)=ky;
tt(kj)=1;



end;
end;
kj=kj+1;
t(j)=1./NM(1,1);
r(j)=NM(2,1)./NM(1,1);
T=(cos(thetaout)./cos(thetain))*nout*abs(t).^2./nin;
R=abs(r).^2;
end;

%% plot results
plot(w,T,'LineWidth',2)
hold on
plot(w,R,'r')
plot(w,T+R,'k--')
xlabel('angular frequency (rad/s)');
ylabel('T(blue),R(red),T+R(black)');
axis([0 maxw 0 1.1])