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

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;