pbgband1.m

本程序可计算二维光子晶体的带隙

function Bandwidth=PBGBand1(ea,eb,f,PCType,Keach,TEorTM,NEIG)

%function PBGBand(ea,eb,R,PCType,Keach,TEorTM)
%--------------------------------------------------------------
%| This is a program to calculate the Photonic Bands of two   |
%| dimension Photonic Crystal with circular inclusions.       |
%| It calculates both TE and TM modes (E and H polarization)  |
%--------------------------------------------------------------
%Parameters:
%ea: The dielectric constant of the circular inclusions.
%eb: The dielectric constant of background.
%R: The radius of dielectric columns
%PCType  =1: Square lattice
%        =2: Triangular lattice
%        =3: Honeycomb
%Keach: The number of k vectors in each wave vector branch.
%TEorTM: =0: TE modes
%        =1: TM modes
%--------------------------------------------------------------



%disp('--------------------------------------------------')
%if (TEorTM==0)
   % disp('Plane wave expansion method for PC bands: TE modes');
   %else
   % disp('Plane wave expansion method for PC bands: TM modes');
   %end    
%disp('--------------------------------------------------')
%if (PCType==1)
    %disp('Square lattice');
    %end
%if (PCType==2)
 %   disp('Triangular lattice');
 %end%
%if (PCType==3)
 %   disp('Honeycomb lattice');
 %end
ea=1;
PCType=1;
Keach=20;
TEorTM=0;
NEIG=10;
eb=6;
R=0.2;

%Control parameters
Ktype=3; % The number of band parts, such as X->M, T->X, ...
NumberK=Ktype*Keach; %The total number of K vector;

%Initial parameters
a=1;   %Lattice constance.
if PCType==1
    a1=a*[1,0];
    a2=a*[0,1];
end

if PCType==2
    a1=a*[1,0];
    a2=a*[0.5,sqrt(3)/2];
end

if PCType==3
    a1=a*[sqrt(3),0];
    a2=a*[sqrt(3)/2,1.5];
end
%a1,a2 are the basic vectors of lacctice cell.

ac=abs(a1(1)*a2(2)-a1(2)*a2(1));
if PCType==1|PCType==2
    f=pi*R*R/ac;
%     R=sqrt(f*ac/pi);
else
    f=2*pi*R*R/ac;
%     R=sqrt(f*ac/(2*pi));
end
%ac: Area of lattice cell.

b1=2*pi/ac*[a2(2),-a2(1)];
b2=2*pi/ac*[-a1(2),a1(1)];
%b1, b2 are vectors in reciprocal space.

%f=pi*R*R/ac;
%f: The filling fraction, i.e. the fraction of
%   the total volume occupied by the rods.

MaxDimForG=9;  % The max Potive Number of the reciprocal lattice, G
DimForG=2*MaxDimForG+1; 
NPW=DimForG*DimForG; %NPW: The number of Plane Waves

%disp('--------------------------------------------------')
%disp('Dielectric constant FT--- BEGIN')

gtemp=-MaxDimForG:MaxDimForG;
gtemp1=repmat(gtemp,DimForG,1);
Gx=b1(1)*gtemp1+b2(1)*gtemp1';
Gy=b1(2)*gtemp1+b2(2)*gtemp1';
Gx=Gx(:)';
Gy=Gy(:)';

%disp(strcat('The number of plane waves is--',num2str(NPW)));

Gx_m=repmat(Gx,NPW,1);
Gx_n=Gx_m';
Gy_m=repmat(Gy,NPW,1);
Gy_n=Gy_m'; 

%calculate the Matrix coefficience.
ek0=f/ea+(1-f)/eb; 
ekc=(1/ea-1/eb)*f*2;  

%Calculate the ek matrix, the coefficence of Fourier transform of ek.
GR_mat=sqrt((Gx_m-Gx_n).*(Gx_m-Gx_n)+(Gy_m-Gy_n).*(Gy_m-Gy_n))*R;
if PCType==1|PCType==2
   %eliminate the division on zero in the calculatation of ek
   na=find(GR_mat==0);
   GR_mat(na)=1;   
   ek_mat=ekc*besselj(1,GR_mat)./GR_mat; 
   ek_mat(na)=ek0;
end
if PCType==3
   %eliminate the division on zero in the calculatation of ek
   na=find(GR_mat==0);
   GR_mat(na)=1;   
   ek_mat=cos((Gx_m-Gx_n).*a/2+(Gy_m-Gy_n).*a*sqrt(3)/6).*ekc.*besselj(1,GR_mat)./GR_mat; 
   ek_mat(na)=ek0;
end
%toc
%tic
%Calculated points:
Point=zeros(Ktype+1,2);
if PCType==1 %Square lattice
   Point(1,:)=[(b1(1)+b2(1))/2,(b1(2)+b2(2))/2]; %M point
   Point(2,:)=[0,0]; %Gama Point
   Point(3,:)=[b1(1)/2,0]; %X Point
   Point(4,:)=[(b1(1)+b2(1))/2,(b1(2)+b2(2))/2]; %M point
end
if PCType==2|PCType==3 %Triangular lattice and Honeycomb lattice
    Point(1,:)=[0,b2(2)/2]; % M point
    Point(2,:)=[0,0]; %Gama Point
    Point(3,:)=[b2(2)*sqrt(3)/6,b2(2)/2]; %K Point
    Point(4,:)=[0,b2(2)/2]; % M point
   %Point(1,:)=[0,0]; %Gama Point
   %Point(2,:)=[0,b2(2)/2]; % M point
   %Point(3,:)=[b2(2)*sqrt(3)/6,b2(2)/2]; %K Point
   %Point(4,:)=[0,0]; %Gama Point
end

%These three are for the K vectors, for the different case.
K1=[];
K2=[];

for ktnum=1:Ktype
   K1temp=linspace(Point(ktnum,1),Point(ktnum+1,1),Keach+1);
   K2temp=linspace(Point(ktnum,2),Point(ktnum+1,2),Keach+1);
   K1=[K1,K1temp(1:Keach)];
   K2=[K2,K2temp(1:Keach)];
end
%disp('Dielectric constant FT--- END')
%disp('--------------------------------------------------')

%disp('Eigen value calculations--- BEGIN')

eigval=[]; %Initial the eigvalue matrix.

for knum=1:NumberK
   %disp(strcat('---K vector No.',num2str(knum),'---',num2str(NumberK))) 
   
   kx=K1(knum);
   ky=K2(knum);
   %tic
   
   %Now begin to calculate the matrix H:    
   if (TEorTM==0)
       %TE part
       KGmn_mat=(kx+Gx_m).*(kx+Gx_n)+(ky+Gy_m).*(ky+Gy_n);
       %KGmn_mat=(kx-Gx_m).*(kx-Gx_n)+(ky-Gy_m).*(ky-Gy_n);
       H=KGmn_mat.*ek_mat;
   else
       KGmn_mat=(kx+Gx_m).^2+(ky+Gy_m).^2;
       %KGmn_mat=(kx-Gx_m).^2+(ky-Gy_m).^2;
       H=KGmn_mat.*ek_mat;
       %TM part
       %Place your codes below for Task 1.   
   end       
       
   %Find the eigenvalues
   eigvalue=sort(eig(H));
   eigval=[eigval,eigvalue(1:NEIG)] ;
end

eigval=[eigval,eigval(:,1)];  
eigval=sqrt(eigval)*a*0.5/pi;
eigval=real(eigval); %Complete All the things
wdiff=[];
for k=1:NEIG-1
    wmax=max(eigval(k,: ),[],2);
    wmin=min(eigval(k+1,:),[],2);
    if (wmin-wmax)>0
        wdiff=[wmax;wmin];
        break
    end
end
if size(wdiff,2)>0
    Bandwidth=wdiff(2)-wdiff(1);
else
    Bandwidth=[];
   return
end