二维光子晶体带隙绘制程序_平面波展开法(最终版)
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;
