inplane1d.m

plane wave expansion for 2D and 3D photonic crystal band gap calculations

% this program is designed for solving the 1D photonic structure
% Structure is given by Dr. Arnel C Lavarias
% The structure is a periodic multilayer structure with 3 medium
% ---------------------------------------------------------------
% 1D multilayer structure with normal incidence
% ---------------------------------------------------------------
% Shangping Guo, 2000
% Old Dominion University 
clear
warning off
epsa=1;
epsb=1.46^2;
epsc=12;
a=0.75; %um
ra=0.4; %um
rb=0.05;%um
rc=0.25/2;%um

a1=a;
b1=2*pi/a;

n=input('please input n: ');
%n=24;
NumberofPW=(2*n+1);

count=1;
G=(-n:n).*b1;
r=(-n:n);
N=4*n+1;
%N=2*NumberofPW+1;
mc=round(N*rc/a);
mb=round(N*rb/a);
ma=N-2*mc-2*mb;
%this is very important to expand the epsilon function like this
%even expansion, so that the imag part is omitted
eps1=[epsc*ones(mc,1);epsb*ones(mb,1);
   epsa*ones(ma,1);epsb*ones(mb,1);epsc*ones(mc,1)];
eps20=real(fftshift(fft(eps1)./N));
%%Now we get the FFT matrix for the dielectric function. We need to get the E(G,G') matrix
%%this matrix include the frequency all we need

for x=1:NumberofPW,
  	for y=x:NumberofPW,
        b=r(x)-r(y)+2*n+1;
        eps2(x,y)=eps20(b);
        eps2(y,x)=eps2(x,y);   
   end
end

k0=-pi/a:2*pi/a/30:pi/a;  %%only kx is considered.
counter=1;
for ii=1:length(k0),
   k=k0(ii);
   M=abs((k+G')*(k+G)).*inv(eps2);
	E=sort(abs(eig(M)));
	freq(:,counter)=sqrt(abs(E(1:10))).*a./2./pi;
	display(sprintf('calculation of k=%f is finished',k));
	counter=counter+1;
end
tmpx=1:length(k0);
plot(tmpx,freq,'linewidth',2)
title('Band structure of a 1D photonic band structure (Arnel)')
xlabel('wave vector')
ylabel('wa/2\pic')
grid on
axis([1 31 0 1])