gapkytm1d.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 off-plane incidence: TE wave
% ---------------------------------------------------------------
% Shangping Guo, 2000
% Old Dominion University 
clear
warning off
epsa=1.4^2;
epsb=3.4^2;
a=1.0; %um
ra=0.5; %um
rb=0.5;%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;

mb=round(N*rb/a/2);
ma=N-2*mb;
%this is very important to expand the epsilon function like this
%even expansion, so that the imag part is omitted
eps1=[epsb*ones(mb,1);epsa*ones(ma,1);epsb*ones(mb,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
i=sqrt(-1);
eps2=inv(eps2);
bandcount=1;
for ky=(0:0.1:1)*2*pi/a,
k0=(-pi/a:2*pi/a/30:pi/a)+i*ky;  %%only kx is considered.
counter=1;
for ii=1:length(k0),
   k=k0(ii);
   %   M=abs((k+G')*(k+G)).*inv(eps2); %in-plane
	M=abs(k+G.')*abs(k+G).*eps2; %TM wave
	%M=(real(k+G.')*real(k+G)+imag(k+G.')*imag(k+G)).*eps2; %%TE polarization
	E=sort(abs(eig(M)));
	freq(:,counter)=sqrt(abs(E)).*a./2./pi;
	display(sprintf('calculation of k=%f is finished',k));
	counter=counter+1;
end
%%%find band gap in these calculated levels and save them
display(sprintf('Searching bandgap for Ky=%f',ky));
for ii=1:5,
   %%find a band gap, save these values
   if ii==1 & min(freq(ii,:)) >= 0
         band(1,bandcount)=ky*a/pi/2;
         band(2,bandcount)=0;%lower level of the band gap
         band(3,bandcount)=min(freq(ii,:));%upper level of the band gap
         band(4,bandcount)=0;%the number of the band gap level
         display(sprintf('find a band (%f %f %f):',band(1,bandcount),band(2,bandcount),band(3,bandcount)));
         bandcount=bandcount+1;
   elseif min(freq(ii,:)) > max(freq(ii-1,:))+0.005   
         band(1,bandcount)=ky*a/pi/2;
         band(2,bandcount)=max(freq(ii-1,:));%lower level of the band gap
         band(3,bandcount)=min(freq(ii,:));%upper level of the band gap
         band(4,bandcount)=ii-1;%the number of the band gap level
         display(sprintf('find a band (%f %f %f):',band(1,bandcount),band(2,bandcount),band(3,bandcount)));
         bandcount=bandcount+1;
   end
end   %%end search for band gap
   
end	%%end the ky for loop   

ab1=band;
count1=1;
while length(band)>0,
   ab(:,2,count1)=band(:,1);
   ab(:,1,count1)=1;
   band(:,1)=[];
   count2=2;
   ii=1;
   while ii<=size(band,2),
      if band(4,ii)==ab(4,2,count1),
         count2=count2+1;
         ab(:,count2,count1)=band(:,ii);
         ab(:,1,count1)=count2-1;
         band(:,ii)=[];
      else
         ii=ii+1;
      end
   end
   count1=count1+1;
end
for ii=1:5
   length=ab(1,1,ii)+1;
   plot(ab(1,2:length,ii),ab(2,2:length,ii),'o-',ab(1,2:length,ii),ab(3,2:length,ii),'-+')
   hold on
end
title('TE polarization with a ky')
xlabel('Normalized ky*a/2\pi')
ylabel('Normalized f=a/\lambda')
hold off
axis([0 1 0 0.8])

%plot(ab1(1,:),ab1(2,:),'-o',ab1(1,:),ab1(3,:),'-+')