honeymap.m

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

% this program is designed for solving the two dimensional photonic structure
% we are using the parameters from one of the paper
% The structure is a triangular 2d periodic structure
clear
warning off
epsa=11.4;
epsb=1;
a=1.0;
%f=0.8358;
%R=sqrt(sqrt(3).*f/(2*pi)*a^2);
i=sqrt(-1);
f=0.5;
R=sqrt(f/(4*pi*sqrt(3)/9))*a;
bandcount=1;
for R=0.2:0.02:0.2,
  
f=4*pi*sqrt(3)/9*(R/a)^2;

%f=2*pi/sqrt(3)*R^2/a^2;
%b1=2*pi/a*[1 -1/sqrt(3) 0];
%b2=2*pi/a*[0 2/sqrt(3) 0];
NumberofCell=1;
%b1=4*pi/sqrt(3)/9/a*(-1.5-sqrt(3)/2*i)/NumberofCell;
%b2=4*pi/sqrt(3)/9/a*(-1.5+sqrt(3)/2*i)/NumberofCell;
b1=4*pi*sqrt(3)/9/a*(1.5+sqrt(3)/2*i)/NumberofCell;
b2=4*pi*sqrt(3)/9/a*(1.5-sqrt(3)/2*i)/NumberofCell;
%n=input('please input n: ');
n=5;
NumberofPW=(2*n+1)^2;

count=1;
for x=-n:n,
  	for y=-n:n,
        G(count)=x*b1+y*b2;
     %   r(count)=x*a1+y*a2;
      count=count+1;
   end
end

for x=1:NumberofPW,
  	for y=x+1:NumberofPW,
        eps2(x,y)=(epsa-epsb)*2*f*besselj(1,abs(G(x)-G(y))*R)./(abs(G(x)-G(y))*R)*...
           cos(imag(G(x)-G(y))*a/2);
   	eps2(y,x)=eps2(x,y);   
   end
  	eps2(x,x)=f*epsa+(1-f)*epsb;
end
   
k1=(0:0.1:1.0)*sqrt(3)/2.*i*4*pi*sqrt(3)/9/a;
k2=((0.1:0.1:1.0)./2+sqrt(3)/2*i).*4*pi*sqrt(3)/9/a;
k3=(0.9:-0.1:0).*(1.0/2.0+sqrt(3)/2*i).*4*pi*sqrt(3)/9/a;%-(1/3+1/sqrt(3)*i)*2*pi/a;
k0=[k1 k2 k3];
%k0=k1;
%k0=(1/3+1/sqrt(3)*i).*2*pi/a;
counter=1;
eps2=inv(eps2);
for ii=1:length(k0),
   k=k0(ii);
   M=(real(k+G.')*real(k+G)+imag(k+G.')*imag(k+G)).*(eps2);%%%For TE mode
  % M=abs(k+G.')*abs(k+G).*(eps2);%This is for TM mode
   E=sort(abs(eig(M)));
	freq(:,counter)=sqrt(abs(E(1:20))).*a./2./pi;
%	display(sprintf('calculation of k=%f+%fi is finished',real(k),imag(k)));
	counter=counter+1;
end
   display(sprintf('Searching bandgap for R=%f',R));
   for ii=1:14,
%%this judgement condition is important, for it may bring a lot of very narrow band gaps
        amin=min(freq(ii+1,:));
        amax=max(freq(ii,:));
        center=(amin+amax)/2;
        width=amin-amax;
        if width>center*0.002

%      if min(freq(ii+1,:)) > max(freq(ii,:))   %%find a band gap, save these values
         band(1,bandcount)=R;
         band(2,bandcount)=max(freq(ii,:));%lower level of the band gap
         band(3,bandcount)=min(freq(ii+1,:));%upper level of the band gap
         band(4,bandcount)=ii;%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
ab1=band;
count1=1;
while length(band)>0,
   ab(:,1,count1)=band(:,1);
   band(:,1)=[];
   count2=1;
   ii=1;
   while ii<=size(band,2),
      %display(sprintf('ii=%f band(4,ii)=%f',ii,band(4,ii)));
      if band(4,ii)==ab(4,1,count1),
         count2=count2+1;
         ab(:,count2,count1)=band(:,ii);
         band(:,ii)=[];
                else
                ii=ii+1;
                end
   end
   arraylength(count1)=count2;
   count1=count1+1;
end
%now we are going to plot the gap map
%%%
for ii=1:size(ab,3),
ll=arraylength(ii);
x=ab(1,1:ll,ii);
y1=ab(2,1:ll,ii);
y2=ab(3,1:ll,ii);
fill([x,rot90(x)'],[y1,rot90(y2)'],'b','erasemode','xor')
hold on
end
axis([0 0.5 0 1])
title('Gap map: Honeycom lattice of GaAs columns')
xlabel('Radius:r/a')
ylabel('wa/2\pic')
text(0.05,0.9,'Red:TE blue:TM')