sqdos2d.m

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

%%         2D PBG:Square lattice
%%         Shangping Guo
%%         sqband.m
clear all
warning off;
epsa=8.9;  %material for the atom
epsb=1;    %material of the background
a=1.0;
R=0.20*a;  %radius of the atom
i=sqrt(-1);
f=pi*R^2/a^2;  %volume fraction of atom(s) in the cell
NumberofCell=1; %supercell size=1X1
a1=a;
a2=a*i;
b1=2*pi/a/NumberofCell;
b2=2*pi/a/NumberofCell*i;
n=input('please input n:');
display('Fourier transforming.....');
tic;
NumberofPW=(2*n+1)^2;
mind=(-2*n:2*n)'+2*n+1;
mind=mind(:,ones(1,size(mind)))-2*n-1;
GG=mind'*b1+mind*b2;
%clear mind;
eps21=2*f*(epsa-epsb)*besselj(1,abs(GG).*R)./(abs(GG).*R);
eps21(2*n+1,2*n+1)=epsb+f*(epsa-epsb);
%zz=[0,0]*[a1 a2].';%use 1X1 supercell to verify the algorithm
%5X5 super cell
zz=[
-2 -2;-2 -1;-2 0;-2 1;-2 2; 
-1 -2;-1 -1;-1 0;-1 1;-1 2;
 0 -2; 0 -1;      0 1; 0,2; 
 1 -2; 1 -1; 1 0; 1 1; 1 2;
 2,-2; 2,-1; 2 0; 2,1; 2,2]*[a1 a2].';
zz=[0,0]*[a1 a2].';%use 1X1 supercell to verify the algorithm
%zz=[0 0; 0 1;1 0;1 1]*[a1 a2].';%2X2 supercell
%zz=[-1 -1;-1 0;-1 1;0 -1;0 0;0 1;1 -1;1 0;1 1]*[a1 a2].';
eps20=zeros(length(eps21));
for x=1:length(zz),
   eps20=eps20+exp(i*(real(GG).*real(zz(x))+imag(GG).*imag(zz(x)))).*eps21;
end
ff=pi*R^2*length(zz)/(NumberofCell^2*a^2);
eps20=eps20./NumberofCell^2;
eps20(2*n+1,2*n+1)=epsb+ff*(epsa-epsb);
count=1;
for y=-n:n,
 for x=-n:n;
      G(count)=x*b1+y*b2;
      r(count,:)=[x,y];
      count=count+1;
   end
end

display('Building eps(G,G) matrix from the FFT matrix');
for x=1:NumberofPW,
  	for y=x:NumberofPW,
	b=r(x,:)-r(y,:)+2*n+1;        
	eps2(x,y)=eps20(b(2),b(1));
   eps2(y,x)=eps2(x,y);   
   end
end
  
  %form the k0 matrix
  mesh=5; %divided the 1Bz into 21X21 mesh
  b10=b1/mesh/2;
  b20=b2/mesh/2;
  k0count=1;
  for ii=-mesh:mesh,
     for jj=-mesh:mesh,
        k0(k0count)=ii*b10+jj*b20;
        k0count=k0count+1;
     end
  end
  
display('Now calculate the eigen values..');
eps2=inv(eps2);
if max(max(real(eps2))) > 10^7*max(max(imag(eps2)))
display('Your lattice is inversion symmetric');
eps2=real(eps2);
else
display('Imaginary part of FFT is not zero');       
stop;
%here we only demonstrate the inversion symmetric case
%%however, nonsymmetric case is also supported
end   

counter=1;
for ii=1:length(k0),
	k=k0(ii);
	%k=k0;
   %M=(real(k+G.')*real(k+G)+imag(k+G.')*imag(k+G)).*eps2; %TE wave
   M=abs(k+G.')*abs(k+G).*eps2; %TM wave
	E=sort(abs(eig(M)));
	freq(:,counter)=sqrt(abs(E(1:20))).*a./2./pi;
	display(sprintf('%d/%d is finished',ii,length(k0)));
	counter=counter+1;
end
toc;
[max(freq(1,:)),min(freq(2,:))]
tmpx=1:length(k0);
plot(tmpx,freq,'o-','linewidth',2)
title('TM:Band structure of a 2D square PBG with a point defect (5X5)')
xlabel('wave vector')
ylabel('wa/2\pic')
grid on
%save sqdos freq;
%now sorting to get the DOS in (0-1)
%actually, the band is partly sorted, we can do the numbering
%band by band, and dividing the frequency 0-1 into 100 sections
step=0.005;
ndos=zeros(max(freq(15,:))/step+10,1);
for band=1:15,
   index=round(freq(band,:)/step+0.5);
   %nods(index)=ndos(index)+1; 
   %doesnot work for some reason, 
   %maybe because of parallel algorithm, have to do it in a loop
   for ii=1:length(index),
      ndos(index(ii))=ndos(index(ii))+1;
   end
   
end
tmpx1=0:step:1;
plot(tmpx1,ndos(1:length(tmpx1))./max(ndos(1:length(tmpx1))),'.-')
%title('2D square lattice: TM band structure using 40401 k-vectors in 1BZ') 
%xlabel('In-plane k vector')
%ylabel('Normalized f:a/\lambda')
grid