matlab仿真二维光子晶体最简程序

本程序为初学者使用，只考虑MT方向

下面的程序为matlab代码

只考虑MT方向

%This is a simple demo for Photonic Crystals simulation 

%This is a simple demo for Photonic Crystals simulation

%This demo is for TE wave only, so only h wave is considered.

%for TM direction only,10 points is considered.

%---------------------------------------M

%| / |

%| / |

%| / |

%| --------------------|X

%| T |

%| |

%| |

%---------------------------------------

%equation :sum_{G',k}(K+G)(K+G')f(G-G')hz(k+G')=(omega/c)^2*hz(k+G)

%G' can considerd as the index of column, and G as index of rows

%[(K+G1)(K+G1)f(G1-G1) (K+G1)(K+G2)f(G1-G2) ][hz(G1)]=(omega/c)^2[hz(G1)]

%[(K+G2)(K+G1)f(G2-G1) (K+G2)(K+G2)f(G2-G2) ][hz(G2)] [hz(G2)]

%or: THETA_TE*Hz=(omega/c)^2*Hz

%by Gao Haikuo

%date:20170411

clear; clc; 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];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%定义晶格的参数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

epsa = 1; %介质柱的介电常数

epsb = 13; %背景的介电常数

Pf = 0.7; %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

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;

T=(2*pi/a)*[epssys 0];

M=(2*pi/a)*[1/2 1/2]; %????????

X=(2*pi/a)*[1/2 0];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%对于简约布里渊区边界上的每个k，求解其特征频率

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

THETA_TE=zeros(NG,NG); %待解的TE波矩阵

Nkpoints=10; %每个方向上取的点数，

stepsize=0:1/(Nkpoints-1):1; %每个方向上的步长

TX_TE_eig = zeros(Nkpoints,NG); %沿MT方向的TE波的待解的特征频率矩阵

XM_TE_eig = zeros(Nkpoints,NG); %沿MT方向的TE波的待解的特征频率矩阵

MT_TE_eig = zeros(Nkpoints,NG); %沿MT方向的TE波的待解的特征频率矩阵

for n=1:Nkpoints %scan the 10 points along the TM direction

fprintf(['\n k-point:',int2str(n),'of',int2str(Nkpoints),'.\n']);

MT_step = stepsize(n)*(T-M)+M; % get the k

%先求非对角线上的元素

for i=1:(NG-1) % G

for j=(i+1):NG % G'

kGi = MT_step+G(i,:); %k+G

kGj = MT_step+G(j,:); %K+G'

THETA_TE(i,j)=f(i,j)*dot(kGi,kGj); %(K+G)(K+G')f(G-G')

THETA_TE(j,i)=conj(THETA_TE(i,j));

end

end

%再求对角线上的元素

for i=1:NG

kGi = MT_step+G(i,:);

THETA_TE(i,i)=f(i,i)*norm(kGi)*norm(kGi);

end

MT_TE_eig(n,:)=sort(sqrt(eig(THETA_TE))).';

end

%draw

kaxis = 0;

TXaxis = kaxis:norm(T-X)/(Nkpoints-1):(kaxis+norm(T-X));

kaxis = kaxis + norm(T-X);

XMaxis = kaxis:norm(X-M)/(Nkpoints-1):(kaxis+norm(X-M));

kaxis = kaxis + norm(X-M);

MTaxis = kaxis:norm(M-T)/(Nkpoints-1):(kaxis+norm(M-T));

kaxis = kaxis + norm(M-T);

Ntraject = 3;

figure (1)

hold on;

Nk=Nkpoints;

for k=1:NG

for i=1:Nkpoints

EigFreq_TE(i+0*Nk) = TX_TE_eig(i,k)/(2*pi/a);

EigFreq_TE(i+1*Nk) = XM_TE_eig(i,k)/(2*pi/a);

EigFreq_TE(i+2*Nk) = MT_TE_eig(i,k)/(2*pi/a);

end

plot(TXaxis(1:Nk),EigFreq_TE(1+0*Nk:1*Nk),'r',...

XMaxis(1:Nk),EigFreq_TE(1+1*Nk:2*Nk),'r',...

MTaxis(1:Nk),EigFreq_TE(1+2*Nk:3*Nk),'r');

end

grid on;

xlabel('K-Space');

yLabel('Frequency(\omegaa/2\piC)');

axis([0 MTaxis(Nkpoints) 0 1]);

set(gca,'XTick',[TXaxis(1), TXaxis(Nkpoints),...

XMaxis(Nkpoints),MTaxis(Nkpoints)]);

xtixlabel = strvcat('T','X','M','T');

set(gca,'XTickLabel',xtixlabel);