本程序为二维光子晶体Matlab仿真程序，该结果与文献【1】Molding the flow of light，p68 figure 2相互吻合
主程序
%This is a simple demo for Photonic Crystals simulation
%10 points is considered.
%---------------------------------------M
%| / |
%| / |
%| / |
%| --------------------|X
%| T |
%| |
%| |
%---------------------------------------
%by Gao Haikuo
%date:20170411
clear; clc; epssys=1.0e-6; %设定一个最小量，避免系统截断误差或除0错误
global NG G f Nkpoints
%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; %每个方向上取的点数，
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%定义晶格的参数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
epsa = 8.9; %inner
epsb = 1; %outer
Pf = 0.1257; %Pf = Ac/Au 填充率，可根据需要自行设定
Au =a^2; %二维格子原胞面积
Rc = (Pf *Au/pi)^(1/2); %介质柱截面半径
Ac = pi*(Rc)^2; %介质柱横截面积
%construct the G list
NrSquare = 10;
NG =(2*NrSquare+1)^2; % NG is the number of the G value
G = zeros(NG,2);
i = 1;
for l = -NrSquare:NrSquare
for m = -NrSquare:NrSquare
