rectan_circl_pbg.m

利用超晶胞和平面波方法计算了长方格子圆柱光子晶体的能带结构

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%     长方格子圆柱能带，
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
opengl neverselect
warning off
tic 
epsa=12;
epsb=1.3;
Ly=1.0;
Lx=1;
R=0.25*Ly;
i=sqrt(-1);
f=pi*R*R/Ly/Lx;
b1=2*pi/Lx;
b2=2*pi/Ly*i;
band=15;
n=5;
Numberofpw=(2*n+1)^2;
mind=(-Numberofpw:Numberofpw)'+Numberofpw+1;
mind=mind(:,ones(1,size(mind)))-Numberofpw-1;
    
    GG=mind'*b1+mind*b2;     %  超晶胞的倒格矢  % 这里的 GG 取得的 G 远大于 (2n+1)^2
   
  eps21=2*f*(1/epsa-1/epsb)*besselj(1,abs(GG)*R)./(abs(GG)*R);
  eps21(Numberofpw+1,Numberofpw+1)=1/epsb+f*(1/epsa-1/epsb);   
    
 count=1;
%r=ones((2*n+1)^2,2);%?????????????????????????????????????????????????
r=ones((2*n+1)^2,2);%某个 G 的h1和h2的值
for y=-n:n
    for x=-n:n
        G(count)=x*b1+y*b2;  %count是指某个 G 的标号
        r(count,:)=[x,y]; % 一个 r 就是一个 h1,h2 的组合,x , y即是 h1和 h2, 一个  r 对应着一个 G，这里的 G 对应于 本征方程中的 全部求和的 G 即上面 GG 中的一部分
        count=count+1;
    end
end
for x=1:Numberofpw
    for y=x:Numberofpw
        bb=r(x,:)-r(y,:)+(2*n+1)^2+1;
        eps2(x,y)=eps21(bb(1),bb(2));
        eps2(y,x)=eps2(x,y);
     end
 end

 %switch  Koftype
 %case  1
k1=2*pi/Lx*0.5*(0:0.1:1);    %T->X
k2=2*pi/Lx*0.5+2*pi/Ly*(0:0.1:1).*0.5*i;%  X->M
k3=2*pi/Lx*0.5*(0.9:-0.1:0)+2*pi/Ly*0.5*i.*(0.9:-0.1:0);  % M->T

k0=[k1 k2 k3 ];
%case 2
    
 k11=2*pi/Lx*0.5*(1:-0.1:0)+2*pi/Ly*0.5*i.*(1:-0.1:0);  % M->T
   
k22=2*pi/Lx*0.5*(0.1:0.1:1);    %T->X
k33=2*pi/Lx*0.5+2*pi/Ly*(0.1:0.1:1).*0.5*i;%  X->M
k44=2*pi/Lx*0.5*(0.9:-0.1:0)+2*pi/Ly*0.5*i;    %M->X'
k55=2*pi/Ly*(0.9:-0.1:0).*0.5*i;%  X'->T

k00=[k11 k22 k33 k44 k55 ];

%end
kkk=k0;
%eps21=inv(eps21);
counter=1;
%freq(1:Numberofpw,1:length(k0))=0;
freq(1:(2*n+1)^2,1:length(kkk))=0;
for ii=1:length(kkk),
    k=kkk(ii);
    M=abs(k+G.')*abs(k+G).*eps2;
    E=sort(abs(eig(M)));
    freq(:,counter)=sqrt(abs(E))*Ly/2/pi;
    counter=counter+1;
 end
tmpx=1:length(kkk);
figure
plot(tmpx,freq(1:band,:),'linewidth',2)
%title('TM:Band structure of a 2D square PBG ','epsa=',epsa,'epsb=',epsb,'the number of plane wave is ',Numberofpw)
titletext=strcat('TM:Band structure of a 2D ellipes PBG (epsa=',num2str(epsa),', epsb=',num2str(epsb),',Numberofpw=',num2str(Numberofpw),')');
xlabel('wave vector')
ylabel('wa/2\pic')
grid on 
h=title(titletext);
set(h,'FontSize',14);
set(gca,'xtick',[]);

axis([1 length(kkk) 0 1.80])
figure
tmp=ifft2(ifftshift(eps21));
pcolor(fftshift(real(tmp))),shading interp,axis image,colorbar