hontmfld.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);
i=sqrt(-1);
f=0.5;
R=sqrt(f/(4*pi*sqrt(3)/9))*a
%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;
k0=sqrt(3)/2*i*4*pi*sqrt(3)/9/a; %at M point
counter=1;

for ii=1:length(k0),
   k=k0(ii);
%   M=(real(k+G.')*real(k+G)+imag(k+G.')*imag(k+G)).*inv(eps2);
	M=abs(k+G.')*abs(k+G).*inv(eps2);
   [V,F]=eig(M);
   E=diag(F);
	%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
[lamda,index]=sort(E);
hg0=V(:,index(1))';
hg1=V(:,index(2))';
for x=(-2:0.02:2).*a,
for y=(-2:0.02:2).*a,
%r=x+i*y;
fld0(x/0.02+101,y/0.02+101)=exp(i*(real(k+G)*x+imag(k+G)*y))*(hg0.*abs(k+G))';
fld1(x/0.02+101,y/0.02+101)=exp(i*(real(k+G)*x+imag(k+G)*y))*(hg1.*abs(k+G))';
%fld1(x*100+151,y*100+151)=exp(i*(real(k+G)*x+imag(k+G)*y))*(hg1.*abs(k+G))';
end
end
subplot(1,2,1),surf(rot90(real(fld0))),shading interp,view(2),axis square,axis off
axis([1 201 1 201])
title('Lower band:D field at M')
subplot(1,2,2),surf(rot90(imag(fld1))),shading interp,view(2),axis square,axis off
axis([1 201 1 201])
title('Upper band:D field at M')