PhotonicCrystal1D.m

% see https://ocw.mit.edu/courses/materials-science-and-engineering/3-024-electronic-optical-and-magnetic-properties-of-materials-spring-2013/lecture-notes/MIT3_024S13_2012lec23.pdf
clear all; close all; clc;
cols={[0.1216 0.2275 0.5765]} ; % a nice line color
 
% setup 
h1 = 2;  % fractional thickness 1 
h2 = 1;  % fractional thickness 2
n1 = 1;  % = sqrt(eps1)
n2 = 3;  % = sqrt(eps1)
Nk0 = 1000; Nkz = 1000; 
k0 = linspace(0,8,Nk0);
kz = linspace(-1,1,Nkz)*pi/(h1+h2);
 
% definitions
k1 = n1*k0;
k2 = n2*k0;
 
% dispersion equation (lengthy, general form)
M11 = exp(-1i*k1*h1).*(cos(k2*h2)-0.5*1i*(k1./k2+k2./k1).*sin(k2*h2));
M22 = exp(1i*k1*h1).*(cos(k2*h2)+0.5*1i*(k1./k2+k2./k1).*sin(k2*h2));
f = bsxfun(@minus, cos(kz*(h1+h2)), (M11 +M22).'/2);
 
% analytical solution
kza = real(acos(cos(k1*h1).*cos(k2*h2) - (n1^2+n2^2)/(2*n1*n2)*sin(k1*h1).*sin(k2*h2))); 
% nasty fix for visualization purposes - only care about real part above zero
kza(abs(kza)<1e-9) = NaN;               
 
 
% --- figure ---
set_figsize(1,10,17)
% visualize the "dispersion equation" magnitude - black where solutions exist
pcolor(kz*(h1+h2),k0,-log(abs(f))); shading interp
caxis([0.01,3])
colormap(flipud(gray))
hold on 
plot([0,0],minmax(k0),'-k') % axis line
% analytical solution of K as function of Omega
plot(real(kza),k0,'--','color',cols{1},'linewidth',3); 
plot(-real(kza),k0,'--','color',cols{1},'linewidth',3); 
%tidy up
hold off
box on; set(gca,'Layer','Top')
xlabel('\itk_xa','Fontsize',12)
ylabel('\it\omega\ita/c','Fontsize',12)