triang2dte.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=1;
epsb=13;
a=1.0;
%f=0.8358;
%R=sqrt(sqrt(3).*f/(2*pi)*a^2);
R=0.48*a;
i=sqrt(-1);
f=2*pi/sqrt(3)*R^2/a^2;
a1=a;
a2=a*(1/2+sqrt(3)/2*i);
b1=2*pi/a*(1-sqrt(3)/3*i);
b2=2*pi/a*2*sqrt(3)/3*i;
n=input('please input n: ');
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);
   	eps2(y,x)=eps2(x,y);   
   end
  	eps2(x,x)=f*epsa+(1-f)*epsb;
end
   
k1=(0:0.1:1.0)/sqrt(3).*i*2*pi/a;
k2=((0.1:0.1:1.0)./3+1/sqrt(3)*i).*2.*pi./a;
k3=(0.9:-0.1:0).*(1.0/3.0+1/sqrt(3)*i).*2.*pi./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);
	E=sort(abs(eig(M)));
	freq(:,counter)=sqrt(abs(E(1:10))).*a./2./pi;
	display(sprintf('calculation of k=%f+%fi is finished',real(k),imag(k)));
	counter=counter+1;
end
tmpx=1:length(k0);
%plot(n,freq,'o'),hold on
plot(tmpx,freq,'linewidth',2)
title('TE:Band structure of a 2D triangular photonic band structure')
xlabel('wave vector')
ylabel('wa/2\pic')
grid on
axis([0,length(k0)+1,0,1.4])