二维光子晶体带隙 占空比vs带隙结构

程序如下

%This is a simple demo for Photonic Crystals simulation 
%10 points is considered.
%by Gao Haikuo 
%date:20170411

clear; clc;
global NG G f  Nkpoints eigenValue modeset kCorner 
global epsa epsb epssys a b1 b2
epssys=1.0e-6; %设定一个最小量，避免系统截断误差或除0错误

%this is the lattice vector and the reciprocal lattice vector
a=1; a1=a*[1 0]; a2=a*[0 1]; 
b1=2*pi/a*[1 0];b2=2*pi/a*[0 1];
Nkpoints=10; %每个方向上取的点数，
modeset=2;% 1:'TE' 2 'TM'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%定义晶格的参数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
epsa = 8.9; %inner
epsb = 1; %outer
Pf = 0.1257; %Pf = Ac/Au 填充率，可根据需要自行设定
figure(1);
kCorner=(2*pi/a)*[epssys 0;1/2 0;1/2 1/2]; %T X M 

for Pf=0.05:0.05:0.785
    fprintf(['
 Pf:',num2str(Pf),'to',num2str(0.785),'.
']);
    Au =a^2; %二维格子原胞面积
    Rc = (Pf *Au/pi)^(1/2); %介质柱截面半径
    Ac = pi*(Rc)^2; %介质柱横截面积
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %get gap
    [G,f]=getGAndf(Pf,Rc);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %get gap
    eigenValue=getFrequency(kCorner);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %get gap
    gap=getGap();
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %绘draw band
    %drawBand(gap);
    for i=1:size(gap,1)
      plot([Pf ,Pf ],[gap(i,1),gap(i,2)]);
      hold on
    end
    axis([0 0.8 0 0.8]);
end