objpbg.m

计算光子晶体的能带投影图

% 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 square 2d periodic structure
% the higher index rod is in the background of air
clear
warning off
epsa=8.9;
epsb=1;
a=1.0;
R=0.2*a;
i=sqrt(-1);
f=pi*R^2/a^2;
a1=a;
a2=a*i;
b1=2*pi/a;
b2=2*pi/a*i;
%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,abs(G(x)-G(y))*R)./(abs(G(x)-G(y))*R);
eps2(y,x)=eps2(x,y);
end
eps2(x,x)=f*epsa+(1-f)*epsb;
end
eps2=inv(eps2);
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=k0(ii);
M=abs(k+G.')*abs(k+G).*(eps2);
E=sort(abs(eig(M)));
freq(:,counter)=sqrt(abs(E(1:20))).*a./2./pi; % 第counter个波矢 ky 对应头20个本征频率值
% display(sprintf('calculation of k=%f+%fi is finished',real(k),imag(k)));
counter=counter+1;
end
min1(count1)=min(freq(1,:));%第COUNT1个kx波矢的对应的所有ky的第一个本征频率的最小值
max1(count1)=max(freq(1,:));%第COUNT1个kx波矢的对应的所有ky的第一个本征频率的最大值
min2(count1)=min(freq(2,:));%第COUNT1个kx波矢的对应的所有ky的第二个本征频率的最小值
max2(count1)=max(freq(2,:));%第COUNT1个kx波矢的对应的所有ky的第二个本征频率的最大值
count1=count1+1;
end
tmpx=(0:0.1:1)/2;
%plot(tmpx,min1,tmpx,max1,tmpx,min2,tmpx,max2,'linewidth',2)
%plot(tmpx,min1,tmpx,max1,tmpx,min2,'linewidth',2)

hold on
fill([tmpx,rot90(tmpx)'],[tmpx,rot90(max2)'],'g','erasemode','none')
%fill([tmpx,rot90(tmpx)'],[min2,rot90(max2)'],'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(min2,tmpx),rot90(max2)'],'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(2),'facecolor','g')
set(h2(3),'facecolor','r')
set(h2(4),'facecolor','b')