1.主程序

%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 填充率，可根据需要自行设定

    Au =a^2; %二维格子原胞面积
    Rc = (Pf *Au/pi)^(1/2); %介质柱截面半径
    Ac = pi*(Rc)^2; %介质柱横截面积
    kCorner=(2*pi/a)*[epssys 0;1/2 0;1/2 1/2]; %T X M
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %get gap
    [G,f]=getGAndf(Pf,Rc);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %get gap
    eigenValue=getFrequency(kCorner);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %get gap
    gap=getGap();
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %绘draw band
    drawBand(gap);

本程序又可以分为以下几个步骤：

• 求解G和f
• 求解本征频率
• 求解光子带隙
• 绘图

2.求解G和f

function [G,f]=getGAndf(Pf,Rc)
global NG G f  Nkpoints eigenValue kCorner modeset a
global epsa epsb epssys b1 b2
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
        G(i,:)=l*b1+m*b2;
        i = i+1;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%生成k空间中的f(Gi-Gj)的值，i,j 从1到NG。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f=zeros(NG,NG);
for i=1:NG
    for j=1:NG
        Gij=norm(G(i,:)-G(j,:));
        if (Gij < epssys)
            f(i,j)=(1/epsa)*Pf+(1/epsb)*(1-Pf);
        else
            f(i,j)=(1/epsa-1/epsb)*Pf*2*besselj(1,Gij*Rc)/(Gij*Rc);
        end;
    end;
end; 

3.求解本征频率

function eigenValue=getFrequency(kCorner)
global  NG Nkpoints modeset
n_kCorner=size(kCorner,1);
kCorner_ex=[kCorner;kCorner(1,:)];
eigenValue=zeros(Nkpoints,NG,n_kCorner,2);
for i_mode=modeset
    for i_kCorner=1:n_kCorner
        eigenValue(:,:,i_kCorner,i_mode)=...
            getEigValue(kCorner_ex(i_kCorner,:),...
            kCorner_ex(i_kCorner+1,:),i_mode);
    end
end

该程序调用了求解一个布里渊边界本征频率的子程序

function [eigValue]=getEigValue(k_begin,k_end,mode)
%mode :1 for TE 2 for TM
global NG G f Nkpoints
THETA=zeros(NG,NG); %待解的TE波矩阵
stepsize=0:1/(Nkpoints-1):1; %每个方向上的步长
TX_TE_eig = zeros(Nkpoints,NG);
for n=1:Nkpoints %scan the 10 points along the direction
    fprintf([‘\n k-point:‘,int2str(n),‘of‘,int2str(Nkpoints),‘.\n‘]);
    step = stepsize(n)*(k_end-k_begin)+k_begin;  % get the k
    for i=1:(NG-1)   % G
        for j=(i+1):NG % G‘
            kGi = step+G(i,:); %k+G
            kGj = step+G(j,:); %K+G‘
            switch mode
                case 1  %TE mode
                    THETA(i,j)=f(i,j)*dot(kGi,kGj); %(K+G)(K+G‘)f(G-G‘)
                case 2 %TH mode
                     THETA(i,j)=f(i,j)*norm(kGi)*norm(kGj);
            end
            THETA(j,i)=conj(THETA(i,j));
        end
    end
    for i=1:NG
        kGi = step+G(i,:);
        THETA(i,i)=f(i,i)*norm(kGi)*norm(kGi);
    end
    eigValue(n,:)=sort(sqrt(eig(THETA))).‘;
end

4.求解光子晶体带隙

function gap=getGap()
global  eigenValue kCorner modeset NG a
gap=[];
band=[];
n_kCorner=size(kCorner,1);
kCorner_ex=[kCorner;kCorner(1,:)];
for mode=modeset
    for i_kCorner=1:n_kCorner
        for k=1:NG
            subLine=eigenValue(:,k,i_kCorner,mode)*a/(2*pi);
            L=min(subLine);H=max(subLine);
            band=mergeBand(band,L,H);
        end
    end
end
if size(band,1)>=2
    gap=[band(1:end-1,2),band(2:end,1)];
end

里面调用了区间合并的函数mergeBand

function band=mergeBand(band,L,H)
flag=1;
i_L=0;%最小值所在行数
flag_L_in=0;%是否在之间
i_H=0;
flag_H_in=0;
i=1;
flag_H_scan=1;
if isempty(band)
    band=[L,H];
    return;
end
while flag
    if L<band(i,1)
        i_L=i;flag=0;
        flag_L_in=0;
    elseif L<=band(i,2)
        i_L=i;flag=0;
        flag_L_in=1;
    else
        i=i+1;
        if i>length(band(:,1))
            flag=0;flag_H_scan=0;
            band=[band;L,H];
        end
    end
end
flag=1;
if flag_H_scan
    while flag
        if H<band(i,1)
        i_H=i;flag=0;
        flag_H_in=0;
        elseif H<=band(i,2)
            i_H=i;flag=0;
            flag_H_in=1;
        else
            i=i+1;
            if i>length(band(:,1))
                flag=0;
                i_H=i;
            end
        end
    end
    %merge
    if i_L==i_H
        if L<band(i_L,1) && H>=band(i_H,1)
            band(i_L,1)=L;
        elseif H<band(i_L,1)
            band=[band(1:i_L-1,:);[L,H];band(i_L:end,:)];
        end
    else
         if  i_H>length(band(:,1))
             band(i_L,1)=min([band(i_L,1),L]);
             band(i_L,2)=H;
             if i_L+1<=length(band(:,1))
                band(i_L+1:end,:)=[];
             end
         elseif  H>=band(i_H,1)
             band(i_L,1)=min([band(i_L,1),L]);
             band(i_L,2)=max([band(i_H,2),H]);
             band(i_L+1:i_H,:)=[];

         else
             band(i_L,1)=min([band(i_L,1),L]);
             band(i_L,2)=H;
             if i_L+1<=i_H-1
                 band(i_L+1:i_H-1,:)=[];
             end
         end
    end
end

5.最后是绘制图形的程序

function drawBand(gap)
global NG G f  Nkpoints eigenValue kCorner modeset a
n_kCorner=size(kCorner,1);
figure (1)
hold on;
colorSet=[‘r‘,‘b‘];
for mode=modeset
    for i_kCorner=1:n_kCorner
        for k=1:NG
            h=plot((i_kCorner-1)*(Nkpoints-1) : i_kCorner*(Nkpoints-1),...
                eigenValue(:,k,i_kCorner,mode)*a/(2*pi),colorSet(mode));
            set(h, ‘linesmoothing‘, ‘on‘);
        end
    end
end
grid on;
xlabel(‘K-Space‘);
yLabel(‘Frequency(\omegaa/2\piC)‘);
xmax=n_kCorner*(Nkpoints-1);
axis([0 xmax 0 0.8]);
set(gca,‘XTick‘,0:(Nkpoints-1):xmax);
xtixlabel = strvcat(‘T‘,‘X‘,‘M‘,‘T‘);
set(gca,‘XTickLabel‘,xtixlabel);
%draw gap
for i=1:size(gap,1)
    fill([0, xmax,xmax,0],[gap(i,1),gap(i,1),gap(i,2),gap(i,2)],‘r‘);
end
hold off;