二维光子晶体能带的程序采用平面波传输法.docx
- 文档编号:5510167
- 上传时间:2023-05-08
- 格式:DOCX
- 页数:15
- 大小:17.35KB
二维光子晶体能带的程序采用平面波传输法.docx
《二维光子晶体能带的程序采用平面波传输法.docx》由会员分享,可在线阅读,更多相关《二维光子晶体能带的程序采用平面波传输法.docx(15页珍藏版)》请在冰点文库上搜索。
二维光子晶体能带的程序采用平面波传输法
pwm法计算二维光子晶体能带的另一个程序
已经分享了一个,这是另一个,编写方法不同,精度更高。
程序完整含有注释,能同时显示TE和TM波的能带结构。
略作修改也可分开显示。
已经分享了一个,这是另一个,编写方法不同,精度更高。
程序完整含有注释,能同时显示TE和TM波的能带结构。
略作修改也可分开显示。
<<隐藏
function[output_args]=Untitled2(input_args)
%UNTITLED2Summaryofthisfunctiongoeshere
%Detailedexplanationgoeshere
%二维光子晶体带结构计算
%计算二维正方格子或三角格子,
%介质圆柱(介电常数a)立于背景(介电常数b)之中
clear;clc;tic;epssys=1.0e-6;%设定一个最小量,避免系统截断误差或除0错误
%定义实际的正空间格子基矢
latticetype=0;%定义格子类型,1计算正方格子,0计算三角格子
switchlatticetype
case1
a=1;
a1=a*[100];
a2=a*[010];
a3=a*[001];
case0
a=1;
a1=a*[sqrt(3)/21/20];
a2=a*[-sqrt(3)/21/20];
a3=a*[001];
end;
%定义晶格的参数
epsa=1;%介质柱的介电常数
epsb=10.5;%背景的介电常数
Pf=0.28274;%Pf=Ac/Au填充率,可根据需要自行设定
Au=norm(cross(a1,a2));%二维格子原胞面积
Rc=(Pf*Au/pi)^(1/2);%介质柱截面半径
Ac=pi*(Rc)^2;%介质柱横截面积
%生成倒格基失
ra11=(2*pi/a)*cross(a2,a3)/dot(a1,cross(a2,a3));
ra22=(2*pi/a)*cross(a3,a1)/dot(a1,cross(a2,a3));
ra1=ra11(1:
2);
ra2=ra22(1:
2);
%选定参与运算的倒空间格矢量,即参与运算的平面波数量。
%设定一个l,m的取值范围,变化l,m即可得出参与运算的平面波集合。
NrSquare=10;%选定k空间的尺度,即l,m(倒格矢G=l*b1+m*b2)的取值范围,NrSquare确定后,使用平面波数目可能为(2*NrSquare+1)^2
G=zeros((2*NrSquare+1)^2,2);%初始化可能使用的倒格矢矩阵
i=1;
forl=-NrSquare:
NrSquare
form=-NrSquare:
NrSquare
G(i,:
)=l*ra1+m*ra2;
i=i+1;
end;
end;
NG=i-1;%实际使用的平面波数目
G=G(1:
NG,:
);
%还有另一种选取的原则,即考虑到园对称的话,则尽可能使平面波波矢的集合均匀分布
%在一个球内,根据这样的原则,可先选取一个比较大的NrSquare,同时确定最大
%G的模,变化l,m,当G的模值不超过最大值时就选入参与运算的集合中。
这样上面一段
%代码可以改写成
%NrSquare=20;
%Gmax=0.5*NrSquare*norm(ra1);
%G=zeros((2*NrSquare+1)^2,2);
%i=1;
%forl=-NrSquare:
NrSquare
%form=-NrSquare:
NrSquare
%Glm=l*ra1+m*ra2;
%if(norm(Glm<=Gmax))
%G(i,:
)=Glm;
%i=i+1;
%end
%end;
%end;
%NG=i-1;
%G=G(1:
NG,:
);
%生成k空间中的f(Gi-Gj)的值,i,j从1到NG。
f=zeros(NG,NG);
fori=1:
NG
forj=1:
NG
Gij=norm(G(i,:
)-G(j,:
));
if(Gij 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; %定义简约布里渊区的各高对称点 switchlatticetype case1 T=(2*pi/a)*[epssys0]; M=(2*pi/a)*[1/21/2]; X=(2*pi/a)*[1/20]; case0 T=(2*pi/a)*[epssys0]; M=(2*pi/a)*[sqrt(3)/30]; K=(2*pi/a)*[sqrt(3)/31/3]; end; %对于简约布里渊区边界上的每个k,求解其特征频率 THETA_TM=zeros(NG,NG);%待解的TM波矩阵 THETA_TE=zeros(NG,NG);%待解的TE波矩阵 Nkpoints=10;%每个方向上取的点数, stepsize=0: 1/(Nkpoints-1): 1;%每个方向上的步长 switchlatticetype case1 TX_TM_eig=zeros(Nkpoints,NG);%沿TX方向的TM波的待解的特征频率矩阵 TX_TE_eig=zeros(Nkpoints,NG);%沿TX方向的TE波的待解的特征频率矩阵 XM_TM_eig=zeros(Nkpoints,NG);%沿XM方向的TM波的待解的特征频率矩阵 XM_TE_eig=zeros(Nkpoints,NG);%沿XM方向的TE波的待解的特征频率矩阵 MT_TM_eig=zeros(Nkpoints,NG);%沿MT方向的TM波的待解的特征频率矩阵 MT_TE_eig=zeros(Nkpoints,NG);%沿MT方向的TE波的待解的特征频率矩阵 case0 TM_TM_eig=zeros(Nkpoints,NG);%沿TM方向的TM波的待解的特征频率矩阵 TM_TE_eig=zeros(Nkpoints,NG);%沿TM方向的TE波的待解的特征频率矩阵 MK_TM_eig=zeros(Nkpoints,NG);%沿MK方向的TM波的待解的特征频率矩阵 MK_TE_eig=zeros(Nkpoints,NG);%沿MK方向的TE波的待解的特征频率矩阵 KT_TM_eig=zeros(Nkpoints,NG);%沿KT方向的TM波的待解的特征频率矩阵 KT_TE_eig=zeros(Nkpoints,NG);%沿KT方向的TE波的待解的特征频率矩阵 end; forn=1: Nkpoints fprintf(['\nk-point: ',int2str(n),'of',int2str(Nkpoints),'.\n']); %对于TX(正方格子)或TM(三角格子)方向上的每个k值, %求解其特征频率 switchlatticetype case1 TX_step=stepsize(n)*(X-T)+T; case0 TM_step=stepsize(n)*(M-T)+T; end; %先求非对角线上的元素 fori=1: (NG-1) forj=(i+1): NG switchlatticetype case1 kGi=TX_step+G(i,: ); kGj=TX_step+G(j,: ); case0 kGi=TM_step+G(i,: ); kGj=TM_step+G(j,: ); end; THETA_TM(i,j)=f(i,j)*norm(kGi)*norm(kGj); THETA_TM(j,i)=conj(THETA_TM(i,j)); THETA_TE(i,j)=f(i,j)*dot(kGi,kGj); THETA_TE(j,i)=conj(THETA_TE(i,j)); end; end; %%%%% %再求对角线上的元素 %%%%% fori=1: NG switchlatticetype case1 kGi=TX_step+G(i,: ); case0 kGi=TM_step+G(i,: ); end; THETA_TM(i,i)=f(i,i)*norm(kGi)*norm(kGi); THETA_TE(i,i)=f(i,i)*norm(kGi)*norm(kGi); end; %求解TX(正方格子)或TM(三角格子)方向上的k矩阵的特征频率 switchlatticetype case1 TX_TM_eig(n,: )=sort(sqrt(eig(THETA_TM))).'; TX_TE_eig(n,: )=sort(sqrt(eig(THETA_TE))).'; case0 TM_TM_eig(n,: )=sort(sqrt(eig(THETA_TM))).'; TM_TE_eig(n,: )=sort(sqrt(eig(THETA_TE))).'; end; %对于XM(正方格子)或MK(三角格子)方向上的每个k值, %求解其特征频率 switchlatticetype case1 XM_step=stepsize(n)*(M-X)+X; case0 MK_step=stepsize(n)*(K-M)+M; end; %先求非对角线上的元素 fori=1: (NG-1) forj=(i+1): NG switchlatticetype case1 kGi=XM_step+G(i,: ); kGj=XM_step+G(j,: ); case0 kGi=MK_step+G(i,: ); kGj=MK_step+G(j,: ); end; THETA_TM(i,j)=f(i,j)*norm(kGi)*norm(kGj); THETA_TM(j,i)=conj(THETA_TM(i,j)); THETA_TE(i,j)=f(i,j)*dot(kGi,kGj); THETA_TE(j,i)=conj(THETA_TE(i,j)); end; end; %再求对角线上的元素 fori=1: NG switchlatticetype case1 kGi=XM_step+G(i,: ); case0 kGi=MK_step+G(i,: ); end; THETA_TM(i,i)=f(i,i)*norm(kGi)*norm(kGi); THETA_TE(i,i)=f(i,i)*norm(kGi)*norm(kGi); end; %求解XM(正方格子)或MK(三角格子)方向上的k矩阵的特征频率 switchlatticetype case1 XM_TM_eig(n,: )=sort(sqrt(eig(THETA_TM))).'; XM_TE_eig(n,: )=sort(sqrt(eig(THETA_TE))).'; case0 MK_TM_eig(n,: )=sort(sqrt(eig(THETA_TM))).'; MK_TE_eig(n,: )=sort(sqrt(eig(THETA_TE))).'; end; %对于MT(正方格子)或KT(三角格子)方向上的每个k值, %求解其特征频率 switchlatticetype case1 MT_step=stepsize(n)*(T-M)+M; case0 KT_step=stepsize(n)*(T-K)+K; end; %先求非对角线上的元素 fori=1: (NG-1) forj=(i+1): NG switchlatticetype case1 kGi=MT_step+G(i,: ); kGj=MT_step+G(j,: ); case0 kGi=KT_step+G(i,: ); kGj=KT_step+G(j,: ); end; THETA_TM(i,j)=f(i,j)*norm(kGi)*norm(kGj); THETA_TM(j,i)=conj(THETA_TM(i,j)); THETA_TE(i,j)=f(i,j)*dot(kGi,kGj); THETA_TE(j,i)=conj(THETA_TE(i,j)); end; end; %再求对角线上的元素 fori=1: NG switchlatticetype case1 kGi=MT_step+G(i,: ); case0 kGi=KT_step+G(i,: ); end; THETA_TM(i,i)=f(i,i)*norm(kGi)*norm(kGi); THETA_TE(i,i)=f(i,i)*norm(kGi)*norm(kGi); end; %求解TX(正方格子)或TM(三角格子)方向上的k矩阵的特征频率 switchlatticetype case1 MT_TM_eig(n,: )=sort(sqrt(eig(THETA_TM))).'; MT_TE_eig(n,: )=sort(sqrt(eig(THETA_TE))).'; case0 KT_TM_eig(n,: )=sort(sqrt(eig(THETA_TM))).'; KT_TE_eig(n,: )=sort(sqrt(eig(THETA_TE))).'; end; end%这个是'forn=1: Nkpoints'(115行)的循环结束 fprintf('\nCalculationTime: %dsec',toc) savepbs2D %以下部分是绘制光子晶体光子带图 %首先将特定方向(正方格子: TX,XM,MT;三角格子: TM,MK,KT)离散化 switchlatticetype case1 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); case0 kaxis=0; TMaxis=kaxis: norm(T-M)/(Nkpoints-1): (kaxis+norm(T-M)); kaxis=kaxis+norm(T-M); MKaxis=kaxis: norm(M-K)/(Nkpoints-1): (kaxis+norm(M-K)); kaxis=kaxis+norm(M-K); KTaxis=kaxis: norm(K-T)/(Nkpoints-1): (kaxis+norm(K-T)); kaxis=kaxis+norm(K-T); end; Ntraject=3;%所需绘制的特定方向的数目 EigFreq_TM=zeros(Ntraject*Nkpoints,1); EigFreq_TE=zeros(Ntraject*Nkpoints,1); figure (1) holdon; Nk=Nkpoints; fork=1: NG fori=1: Nkpoints switchlatticetype case1 EigFreq_TM(i+0*Nk)=TX_TM_eig(i,k)/(2*pi/a); EigFreq_TM(i+1*Nk)=XM_TM_eig(i,k)/(2*pi/a); EigFreq_TM(i+2*Nk)=MT_TM_eig(i,k)/(2*pi/a); 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); case0 EigFreq_TM(i+0*Nk)=TM_TM_eig(i,k)/(2*pi/a); EigFreq_TM(i+1*Nk)=MK_TM_eig(i,k)/(2*pi/a); EigFreq_TM(i+2*Nk)=KT_TM_eig(i,k)/(2*pi/a); EigFreq_TE(i+0*Nk)=TM_TE_eig(i,k)/(2*pi/a); EigFreq_TE(i+1*Nk)=MK_TE_eig(i,k)/(2*pi/a); EigFreq_TE(i+2*Nk)=KT_TE_eig(i,k)/(2*pi/a); end; end;%fork=1: Nkpoints的结束 switchlatticetype case1 plot(TXaxis(1: Nk),EigFreq_TM(1+0*Nk: 1*Nk),'b',... XMaxis(1: Nk),EigFreq_TM(1+1*Nk: 2*Nk),'b',... MTaxis(1: Nk),EigFreq_TM(1+2*Nk: 3*Nk),'b'); 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'); case0 plot(TMaxis(1: Nk),EigFreq_TM(1+0*Nk: 1*Nk),'b',... MKaxis(1: Nk),EigFreq_TM(1+1*Nk: 2*Nk),'b',... KTaxis(1: Nk),EigFreq_TM(1+2*Nk: 3*Nk),'b'); plot(TMaxis(1: Nk),EigFreq_TE(1+0*Nk: 1*Nk),'r',... MKaxis(1: Nk),EigFreq_TE(1+1*Nk: 2*Nk),'r',... KTaxis(1: Nk),EigFreq_TE(1+2*Nk: 3*Nk),'r'); end; end;%fork=1: NG的结束 gridon; holdoff; titlestr='2DPhotonicbandstructurefor'; switchlatticetype case1 titlestr=[titlestr,'square']; case0 titlestr=[titlestr,'triangle']; end; titlestr=[titlestr,'latticeofdielectriccolumns.']; titlestr=[titlestr,'(epsa=',num2str(epsa),'epsb=',... num2str(epsb),')']; title(titlestr); xlabel('k-Space'); ylabel('Frequency(\omegaa/2\piC)'); switchlatticetype case1 axis([0MTaxis(Nkpoints)01]); set(gca,'XTick',[TXaxis (1)... TXaxis(Nkpoints)... XMaxis(Nkpoints)... MTaxis(Nkpoints)]); xtixlabel=strvcat('T','X','M','T'); case0 axis([0KTaxis(Nkpoints)01]); set(gca,'XTick',[TMaxis (1)... TMaxis(Nkpoints)... MKaxis(Nkpoints)... KTaxis(Nkpoints)]); xtixlabel=strvcat('T','M','K','T'); end; set(gca,'XTickLabel',xtixlabel); end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 二维 光子 晶体 能带 程序 采用 平面波 传输