honeytm.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
%%Honeycomb structure
clear
warning off
epsa=11.4;
epsb=1;
%epsb=1.46^2; %%%Silica glass air system
a=1.0;
%f=0.70;
%R=sqrt(sqrt(3).*f/(2*pi)*a^2);
%R=0.48;
R=0.1;
f=4*pi*sqrt(3)/9*(R/a)^2;
%f=0.5;
%R=sqrt(f/(4*pi*sqrt(3)/9))*a
i=sqrt(-1);
%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,0]/NumberofCell;
b2=4*pi*sqrt(3)/9/a*[1.5,-sqrt(3)/2,0]/NumberofCell;
%n=input('please input n: ');
n=5;

%%%%%%%%%%Adding to test fourier transform
%NumberofPW=(2*n*NumberofCell+1)^2;
%maxPW=2*NumberofPW+1;
%mind=(-NumberofPW:NumberofPW)'+NumberofPW+1;
%mind=mind(:,ones(1,size(mind)))-NumberofPW-1;
%GG=mind*b1+mind'*b2;
%eps21=2*f*(epsa-epsb)*besselj(1,abs(GG).*R)./(abs(GG).*R)*2.*cos(imag(GG)*a/2);
%eps21(NumberofPW+1,NumberofPW+1)=epsb+f*(epsa-epsb);
%clear mind;
%pause
%%%%%%%%%%%%end of adding part
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
%note: in honeycomb structure, we have two rods in one unit cell, so the Fourier transformation
%%of this structure should be modified
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)...
           .*cos(dot(G(x,:)-G(y,:),[0 a/2 0]));
   	eps2(y,x)=eps2(x,y);   
   end
  	eps2(x,x)=f*epsa+(1-f)*epsb;
end
%pause   
%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)=4*pi/a*sqrt(3)/9*sqrt(3)/2.*(0:0.1:1)';
mm0=length(0:0.1:1);
mm1=length(0.1:0.1:1);
k0(mm0+1:mm0+mm1,1)=4*pi*sqrt(3)/a/9*0.5.*(0.1:0.1:1)';
k0(mm0+1:mm0+mm1,2)=4*pi/a*sqrt(3)/9*sqrt(3)/2;
mm0=mm0+mm1+1;
mm1=length(0.9:-0.1:0);
k0(mm0:mm0+mm1-1,1)=4*pi/a*sqrt(3)/9*0.5.*(0.9:-0.1:0)';
k0(mm0:mm0+mm1-1,2)=4*pi/a*sqrt(3)/9*sqrt(3)/2.*(0.9:-0.1:0)';
k0(:,3)=0;%2*pi/a;

%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
   %%%to save computing time, use absk=modulus(K)
   absk=modulus(K);
   e1=[K(:,2)./absk,-K(:,1)./absk,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
   absk=modulus(K);
   e2=cross(e1, K);
   e2=[e2(:,1)./absk,e2(:,2)./absk,e2(:,3)./absk];
   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=([absk.*e2(:,1),absk.*e2(:,2),absk.*e2(:,3)]*[absk.*e2(:,1),absk.*e2(:,2),absk.*e2(:,3)]').*eps2;
   M2=([absk.*e1(:,1),absk.*e1(:,2),absk.*e1(:,3)]*[absk.*e2(:,1),absk.*e2(:,2),absk.*e2(:,3)]').*eps2;
   M3=([absk.*e2(:,1),absk.*e2(:,2),absk.*e2(:,3)]*[absk.*e1(:,1),absk.*e1(:,2),absk.*e1(:,3)]').*eps2;
   M4=([absk.*e1(:,1),absk.*e1(:,2),absk.*e1(:,3)]*[absk.*e1(:,1),absk.*e1(:,2),absk.*e1(:,3)]').*eps2;
%   M=[M1,M3;M2,M4];
	E=sort(abs(eig(M1)));
	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,'linewidth',2)
title('Full band structure of a 2D triangular silica-air hole with a kz')
xlabel('in-plane wave vector k_/_/')
ylabel('wa/2\pic')
grid on
%axis([0,length(k0)+1,0,1.4])