bandkz.m

2d photonic band structure calculation program

% this program is designed for solving the two dimensional photonic structure
% we are using the parameters from one of the paper
clear
warning off
epsa=1;
epsb=13;
%epsb=1.46^2; %%%Silica glass air system
a=1.0;
%f=0.39;
%R=sqrt(sqrt(3).*f/(2*pi)*a^2);
R=0.45;
i=sqrt(-1);
f=2*pi/sqrt(3)*R^2/a^2;
%b1=2*pi/a*(1-sqrt(3)/3*i);
%b2=2*pi/a*2*sqrt(3)/3*i;
b1=2*pi/a*[1 -1/sqrt(3) 0];
b2=2*pi/a*[0 2/sqrt(3) 0];
%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;
      count=count+1;
   end
end

for x=1:NumberofPW,
  	for y=x+1:NumberofPW,
   	eps2(x,y)=(epsa-epsb)*2*f*besselj(1,norm(G(x,:)-G(y,:))*R)./(norm(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;
k0=zeros(length(0:0.1:1)+length(0.1:0.1:1)+length(0.9:-0.1:0),3);
mm0=length(0:0.1:1);
k0(1:mm0,2)=2*pi/a/sqrt(3).*(0:0.1:1)';
mm0=length(0:0.1:1);
mm1=length(0.1:0.1:1);
k0(mm0+1:mm0+mm1,1)=2*pi/a/3.*(0.1:0.1:1)';
k0(mm0+1:mm0+mm1,2)=2*pi/a/sqrt(3);
mm0=mm0+mm1+1;
mm1=length(0.9:-0.1:0);
k0(mm0:mm0+mm1-1,1)=2*pi/a/3.*(0.9:-0.1:0)';
k0(mm0:mm0+mm1-1,2)=2*pi/a/sqrt(3).*(0.9:-0.1:0)';
k0(:,3)=0.0;

%k0=[1/3 1/sqrt(3) 9]*2*pi/a;
counter=1;
eps2=inv(eps2);
for ii=1:size(k0,1),
   k=k0(ii,:);
   K(:,1)=k(1)+G(:,1);
   K(:,2)=k(2)+G(:,2);
   K(:,3)=0;
   %%%find the unit cell perpendicular to it in xy plane
   %%%be sure to deal with the case of modulus(k)=0
   %%%NaN in this case
   e1=[K(:,2)./modulus(K),-K(:,1)./modulus(K),zeros(length(K),1)];
   e1(isnan(e1))=1/sqrt(2); %%%when Kz is not zero,kx,ky is zero choose arbitrary e1
   K(:,3)=k(3)+G(:,3);
   %%%find the other perpendicular unit cell
   %%%be sure to deal with the case of modulus(k)=0
   %%%NaN in this case
   e2=cross(e1, K);
   e2=[e2(:,1)./modulus(e2),e2(:,2)./modulus(e2),e2(:,3)./modulus(e2)];
   e2(isnan(e2))=0;
   %%%form the equation matrix, it should be 2N by 2N
   %%%in this case we will have no TE and TM decoupling
   M1=([modulus(K).*e2(:,1),modulus(K).*e2(:,2),modulus(K).*e2(:,3)]*[modulus(K).*e2(:,1),modulus(K).*e2(:,2),modulus(K).*e2(:,3)]').*eps2;
   M2=([modulus(K).*e1(:,1),modulus(K).*e1(:,2),modulus(K).*e1(:,3)]*[modulus(K).*e2(:,1),modulus(K).*e2(:,2),modulus(K).*e2(:,3)]').*eps2;
   M3=([modulus(K).*e2(:,1),modulus(K).*e2(:,2),modulus(K).*e2(:,3)]*[modulus(K).*e1(:,1),modulus(K).*e1(:,2),modulus(K).*e1(:,3)]').*eps2;
   M4=([modulus(K).*e1(:,1),modulus(K).*e1(:,2),modulus(K).*e1(:,3)]*[modulus(K).*e1(:,1),modulus(K).*e1(:,2),modulus(K).*e1(:,3)]').*eps2;
   M=[M1,M2;M3,M4];

	E=sort(abs(eig(M)));
	freq(:,counter)=sqrt(abs(E(1:20))).*a./2./pi;
	display(sprintf('calculation of k=[%f,%f,%f] is finished',k(1),k(2),k(3)));
	counter=counter+1;
end
tmpx=1:length(k0);
%plot(n,freq,'o'),hold on
plot(tmpx,freq*2*pi,'linewidth',2)
title('Full band structure of a 2D triangular photonic band structure')
xlabel('in-plane wave vector k_/_/')
ylabel('wa/2\pic')
grid on
%axis([0,length(k0)+1,0,1.4])