defectmode.m

利用平面波展开的方法计算二维光子晶体的缺陷模

%%         2D PBG:Square lattice
%%         Shangping Guo
%%         sqband.m
clear all
warning off
tic;
epsa=1;
epsb=1.46.^2;
a=1.0;
f=0.7;
R=sqrt(f/pi*a^2);
%R=0.7*a;;
i=sqrt(-1);
%f=pi*R^2/a^2;
NumberofCell=5;
a1=a;
a2=a*i;
b1=2*pi/a/NumberofCell;
b2=2*pi/a/NumberofCell*i;
%n=input('please input n:');
n=5;
display('Fourier transforming.....');
NumberofPW=(2*n+1)^2;
mind=(-NumberofPW:NumberofPW)'+NumberofPW+1;
mind=mind(:,ones(1,size(mind)))-NumberofPW-1;
GG=mind'*b1+mind*b2;
%clear mind;
eps21=2*f*(epsa-epsb)*besselj(1,abs(GG).*R)./(abs(GG).*R);
eps21(NumberofPW+1,NumberofPW+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].';
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(NumberofPW+1,NumberofPW+1)=epsb+ff*(epsa-epsb);

clear GG;
count=1;
r=ones(2*n*NumberofCell+1,2);
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)^2+1;
        eps2(x,y)=eps20(b(1),b(2));
        eps2(y,x)=eps2(x,y);   
   end
end
  
%k1=2*pi/a*0.5.*(0:0.1:1)*i;;
%k2=2*pi/a*0.5.*(0.1:0.1:1)+2*pi/a*0.5*i;
%k3=2*pi/a*(0.5+0.5*i).*(0.9:-0.1:0);
%kk=[k1 k2 k3]./NumberofCell;
%k0(:,1)=real(kk)';k0(:,2)=imag(kk)';k0(:,3)=0;
GG(:,1)=real(G)';GG(:,2)=imag(G)';GG(:,3)=0;
clear G
G=GG;
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   

%k0(:,3)=13;
k(3)=4.1;

count1=1;
for aa=(0:0.1:1)*pi/a,
bb=(0:0.1:1)*pi/a;
k0=aa+bb*i;

counter=1;
for ii=1:length(k0),
   k(1)=real(k0(ii));k(2)=imag(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(M)));
	freq(:,counter)=sqrt(abs(E)).*a./2./pi;
	display(sprintf('calculation of k=[%f %f %f] is finished',k(1),k(2),k(3)));
	counter=counter+1;
end
	min1(count1)=min(freq(1,:));
	max1(count1)=max(freq(1,:));
	min2(count1)=min(freq(2,:));
   max2(count1)=max(freq(2,:));
   min3(count1)=min(freq(3,:));
	max3(count1)=max(freq(3,:));
   abc(count1,:,:)=freq(:,:);
   count1=count1+1;
end

%tmpx=(0:0.1:1)/2;
%plot(tmpx,min1,tmpx,max1,tmpx,min2,tmpx,max2,tmpx,min3,tmpx,max3,'linewidth',2)
%hold on
%fill([tmpx,rot90(tmpx)'],[tmpx,rot90(max3)'],'g','erasemode','none')
%fill([tmpx,rot90(tmpx)'],[min3,rot90(max3)'],'r','erasemode','none')
%fill([tmpx,rot90(tmpx)'],[min1,rot90(max1)'],'r','erasemode','none')
%fill([tmpx,rot90(tmpx)'],[min(max1,tmpx),rot90(max1)'],'b','erasemode','none')
%fill([tmpx,rot90(tmpx)'],[max(min3,tmpx),rot90(max3)'],'b','erasemode','none')
%axis([0 0.5 0 0.5])
%xlabel('ka/2\pi')
%ylabel('wa/2\pic or a/\lambda')
%title('The projected band for constant-x surface of a square Al rod')
%[h1 h2]=legend('ED states','DE states','EE states')
%set(h2(1),'facecolor','g')
%set(h2(2),'facecolor','r')
%set(h2(3),'facecolor','b')

% i means kx, j means ky
for i=1:11
   for j=1:11
band1(i,j)=abc(i,1,j);
band2(i,j)=abc(i,2,j);
band3(i,j)=abc(i,3,j);
band4(i,j)=abc(i,4,j);
band5(i,j)=abc(i,5,j);
band6(i,j)=abc(i,6,j);
band7(i,j)=abc(i,7,j);
band8(i,j)=abc(i,8,j);
band9(i,j)=abc(i,9,j);
band10(i,j)=abc(i,10,j);
band11(i,j)=abc(i,11,j);
band12(i,j)=abc(i,12,j);
band13(i,j)=abc(i,13,j);
band14(i,j)=abc(i,14,j);
band15(i,j)=abc(i,15,j);
band16(i,j)=abc(i,16,j);
band17(i,j)=abc(i,17,j);
band18(i,j)=abc(i,18,j);
band19(i,j)=abc(i,19,j);
band20(i,j)=abc(i,20,j);
	end
end