平面波法用于一维光子晶体

本算法进一步深入熟悉平面波展开法

结果如下：与文献molding the flow of light 相一致

ps：有多余线，需要进一步求证

主程序如下：

%This is a simple demo for 1D 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 = 13; %inner
epsb = 12; %outer

% Pf = 0.1257; %Pf = Ac/Au 填充率，可根据需要自行设定
% Au =a^2; %二维格子原胞面积
% Rc = (Pf *Au/pi)^(1/2); %介质柱截面半径
% Ac = pi*(Rc)^2; %介质柱横截面积
kCorner=(2*pi/a)*[0.5 0;epssys 0]; %T X M
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%get gap
[G,f]=getGAndf_1D(0.5);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%get gap
gap=getGap();

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%绘draw band
drawBand(gap,0.3);

其中的获取展开系数的getGAndf_1D函数如下

function [G,f]=getGAndf_1D(d)
%这个生成一维空间中的G和f
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=G(i,:)-G(j,:);
        if abs(Gij(2))>epssys
            f(i,j)=0;
        elseif abs(Gij(1))<epssys
            f(i,j)=1/epsb+(1/epsa-1/epsb)*d;
        else
            f(i,j)=(1/epsa-1/epsb)*2/Gij(1)*sin(Gij(1)*d/2);
        end;
    end;
end;