pbgband_2d.m

两个平面波展开程序

% This is a program to calculate the Photonic Bands of two dimension
% Photonic Crystal.
clear

ea = 1;
eb = 13;% the dielectric constant of the background material
R = 0.48;
veinthickness = 0.165;
% Photonic crystals molding the flow of light second edition 
% Chapter 5 figure 5
PCType = 2;
Keach = 6;
TEorTM = 1;
%0 TE modes
%1 TM modes

%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
%--------------------------------------------------------------

tic;

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

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

NEIG=20; %NEIG: The number cut of the Eigvalue.

%Initial parameters
a=1;   %Lattice constance.

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

ac=abs(a1(1)*a2(2)-a1(2)*a2(1)); 
%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.

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

%
%            ^ 
%            | Y
%      O  O  O  O  O  -->(MaxDimForG,MaxDimForG)for this point!
%      O  O  O  O  O 
%    --O--O--O--O--O--> X
%      O  O  O  O  O
%      O  O  O  O  O 
%            |   
%            |
%initial the G matrix
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.
% f=pi*R*R/ac;
% %f: The filling fraction, i.e. the fraction of
% %   the total volume occupied by the rods.
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 4096;
ic = N/2;
jc = N/2;
eps1 = zeros(N,N);
if PCType==1 || PCType==2
    xdist_temp = zeros(1,5);
    ydist_temp = zeros(1,5);
    dist_temp = zeros(1,5);
    for j = 1:N
        for i = 1:N
            xdist_temp(1) = (i - ic)*a1(1) + (j - jc)*a2(1);
            ydist_temp(1) = (i - ic)*a1(2) + (j - jc)*a2(2);
            dist_temp(1) = sqrt((xdist_temp(1)/N)^2 + (ydist_temp(1)/N)^2);
            xdist_temp(2) = (i - ic - N)*a1(1) + (j - jc)*a2(1);
            ydist_temp(2) = (i - ic - N)*a1(2) + (j - jc)*a2(2);
            dist_temp(2) = sqrt((xdist_temp(2)/N)^2 + (ydist_temp(2)/N)^2);
            xdist_temp(3) = (i - ic + N)*a1(1) + (j - jc)*a2(1);
            ydist_temp(3) = (i - ic + N)*a1(2) + (j - jc)*a2(2);
            dist_temp(3) = sqrt((xdist_temp(3)/N)^2 + (ydist_temp(3)/N)^2);
            xdist_temp(4) = (i - ic)*a1(1) + (j - jc - N)*a2(1);
            ydist_temp(4) = (i - ic)*a1(2) + (j - jc - N)*a2(2);
            dist_temp(4) = sqrt((xdist_temp(4)/N)^2 + (ydist_temp(4)/N)^2);
            xdist_temp(5) = (i - ic)*a1(1) + (j - jc + N)*a2(1);
            ydist_temp(5) = (i - ic)*a1(2) + (j - jc + N)*a2(2);
            dist_temp(5) = sqrt((xdist_temp(5)/N)^2 + (ydist_temp(5)/N)^2);
            dist = min(dist_temp);
            if (dist <= R)
                eps1(i,j) = 1/ea;
            else
                eps1(i,j) = 1/eb;
            end
        end
    end
end
if PCType==4
    dist=(1-veinthickness)/2;
    for j=1:N
        for i=1:N
            xdist_temp = abs(i-ic)/N;
            ydist_temp = abs(j-jc)/N;
            if xdist_temp>dist || ydist_temp>dist
                eps1(i,j) = 1/eb;
            else
                eps1(i,j) = 1/ea;
            end
        end
    end
end
%eps1_shift = fftshift(eps1);
eps1_shift = eps1;
eps20 = fftshift(fft2(eps1_shift)/N^2);
ek_mat = zeros(NPW,NPW);
for jj = 1:NPW
    ic_new = ic + floor((jj-1)/DimForG)+1;
    jc_new = jc + rem((jj-1),DimForG)+1;
    for ii = 1:NPW
        ii_new = ic_new - floor((ii-1)/DimForG);
        jj_new = jc_new - rem((ii-1),DimForG);
        ek_mat(ii,jj) = eps20(ii_new,jj_new);
    end
end
ek_mat = real(ek_mat);
%figure,mesh(ek_mat);
%toc
%tic
%Calculated points:
Point=zeros(Ktype+1,2);
if PCType==1 || PCType==4%Square lattice, or rectangular lattice
   Point(1,:)=[0,0]; %Gama Point
   Point(2,:)=[b1(1)/2,0]; %X Point
   Point(3,:)=[(b1(1)+b2(1))/2,(b1(2)+b2(2))/2]; %M point
   Point(4,:)=[0,0]; %Gama Point
end
if PCType==2 || PCType==3 %Triangular lattice and Honeycomb lattice
    Point(1,:)=[0,0]; %Gama Point
    Point(3,:)=[b1(1)/3,b1(2)]; %K Point
    Point(2,:)=[0,b1(2)]; %M 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
       %TM part
       %Place your codes below for Task 1.
       KGmn_mat = sqrt((kx+Gx_m).^2+(ky+Gy_m).^2).*sqrt((kx+Gx_n).^2+(ky+Gy_n).^2);
       H = KGmn_mat.*ek_mat;
   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

max_temp = max(eigval,[],2);
min_temp = min(eigval,[],2);
max_temp1 = [];
min_temp1 = [];
gapnum = 0;
for i = 1:NEIG
    flag = 1;
    for j = 1:NEIG
        if max_temp(i)>min_temp(j) && max_temp(i)<max_temp(j)
            flag = 0;
            break;
        end
    end
    if flag == 1
        gapnum = gapnum + 1;
        max_temp1(gapnum) = max_temp(i);
    end
end
gapnum = 0;
for i = 1:NEIG
    flag = 1;
    for j = 1:NEIG
        if min_temp(i)>min_temp(j) && min_temp(i)<max_temp(j)
            flag = 0;
            break;
        end
    end
    if flag == 1
        gapnum = gapnum + 1;
        min_temp1(gapnum) = min_temp(i);
    end
end
max_temp1 = sort(max_temp1);
min_temp1 = sort(min_temp1);

disp('Eigen value calculation0s--- END')
disp('--------------------------------------------------')

%Plot the figures
%x=linspace(0,10,NumberK+1);
for m=1:Ktype
   D(m)=sqrt((Point(m+1,1)-Point(m,1))^2+(Point(m+1,2)-Point(m,2))^2);
   xtemp(m,:)=linspace(0,D(m),Keach+1);
end
x=xtemp(1,1:Keach);
Dtotal=0;
for m=2:Ktype
   Dtotal=Dtotal+D(m-1);
   x=[x,xtemp(m,1:Keach)+Dtotal];
end
x=[x,xtemp(Ktype,Keach+1)+Dtotal];
x=x/max(x);

MaxB=0.8;
x1=x(Keach+1);
x2=x(Keach*2+1);

figure;
clf;

h=plot(x,eigval,'b-',[x1 x1],[0 MaxB],'k:',[x2 x2],[0 MaxB],'k:');
set(h,'LineWidth',2.0);
for i = 1:length(max_temp1)-1
    X=[0,1,1,0];
    Y=[max_temp1(i),max_temp1(i),min_temp1(i+1),min_temp1(i+1)];
    if min_temp1(i+1)-max_temp1(i) > 0.01
        h=patch(X,Y,[1 1 0]);
    end
end
if (TEorTM==0)
    text(x(Keach*Ktype+1)*0.83,MaxB*0.1,...
        'TE modes',...
        'BackgroundColor',[.7 .9 .7],...
        'EdgeColor','red',...
        'LineWidth',1);
    %legend('TE modes',4);
else
    text(x(Keach*Ktype+1)*0.83,MaxB*0.1,...
        'TM modes',...
        'BackgroundColor',[.7 .9 .7],...
        'EdgeColor','red',...
        'LineWidth',1);
    %legend('TM modes',4);
end

axis([0 1 0 MaxB]);
%h=ylabel('Normalized frequency (a/\lambda)');
ylabel('Normalized frequency (\omegaa/2\pic)');
%set(h,'FontSize',14);
if (PCType==1)
   titletext=strcat('Square Lattice (ea=',num2str(ea),', eb=',num2str(eb),', R=',num2str(R),')');
   %title('Square Lattice (ea=',num2str(ea),', eb=',num2str(eb),', R=',num2str(R),')');
   text(x(1)-0.00,-0.03, '\Gamma')
   text(x1-0.01,-0.03, 'X','FontSize',12)
   text(x2-0.01,-0.03, 'M','FontSize',12)
   text(x(Keach*Ktype+1)-0.00,-0.03, '\Gamma')
end
if (PCType==2)
   titletext=strcat('Triangular Lattice (ea=',num2str(ea),', eb=',num2str(eb),', R=',num2str(R),')');
   text(x(1)-0.00,-0.03, '\Gamma')
   text(x1-0.01,-0.03, 'M','FontSize',12)
   text(x2-0.01,-0.03, 'K','FontSize',12)
   text(x(Keach*Ktype+1)-0.00,-0.03, '\Gamma')
end
if (PCType==3)
   titletext=strcat('Honeycomb Lattice (ea=',num2str(ea),', eb=',num2str(eb),', R=',num2str(R),')');
   text(x(1)-0.02,-0.03, '\Gamma','FontSize',14)
   text(x1-0.02,-0.03, 'K','FontSize',14)
   text(x2-0.02,-0.03, 'M','FontSize',14)
   text(x(Keach*Ktype+1)-0.02,-0.03, '\Gamma','FontSize',14)
end
if (PCType==4)
   titletext=strcat('Square Lattice (ea=',num2str(ea),', eb=',num2str(eb),...
   ', vein thickness =',num2str(veinthickness),')');
   %title('Square Lattice (ea=',num2str(ea),', eb=',num2str(eb),', R=',num2str(R),')');
   text(x(1)-0.00,-0.03, '\Gamma')
   text(x1-0.01,-0.03, 'X','FontSize',12)
   text(x2-0.01,-0.03, 'M','FontSize',12)
   text(x(Keach*Ktype+1)-0.00,-0.03, '\Gamma')
end
h=title(titletext);
%set(h,'FontSize',14);
set(gca,'xtick',[]);

%Save the data
if (TEorTM==0)
    save datate.mat x ea eb R PCType MaxB Keach eigval;
else
    save datatm.mat x ea eb R PCType MaxB Keach eigval;
end    

size(eigval)
toc;