%% 1D photonic band 
clear
close all
%% frequency range
maxw=2.5e15;
minw=1e11;
Nw=2000;
w=linspace(minw,maxw,Nw);
lam0=2*pi*3e14./w;     % vacuum wavelength in micron
kj=1;
for j=1:Nw;
%% unit cell  (AB)%%%%%%%%

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

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

a=dA+dB;               % unit cell thickness

ky=0;
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))
%% 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;
%% homogenous medium
nmed=(nA*dA+nB*dB)./a;  % average index
kxM=2*pi.*nmed./lam0(j);
PNMED=[ exp(-i.*kxM.*a)   0;
      0                  exp(i.*kxM.*a)]; 
%% eigen values %%
Ka=acos(0.5*(M(1,1)+M(2,2)))./a;
KaM(j)=acos(0.5*(PNMED(1,1)+PNMED(2,2)))./a;
    if abs(real(Ka))>0;
    if abs(real(Ka))<pi./a;
K(kj)=real(Ka);
freq(kj)=w(j);
KY(kj)=ky;
K(kj+1)=-Ka;
freq(kj+1)=w(j);
KY(kj+1)=ky;
kj=kj+2;
    end;
    end;

end;

%% plot results
plot(K*a./pi,freq,'.','MarkerSize',4)
hold on
plot(KaM*a./pi,w,'r.','MarkerSize',3)
plot(-KaM*a./pi,w,'r.','MarkerSize',3)

grid on
xlabel('K (units of \pi/a)')
ylabel('angular frequency (rad/s)')
hold off