二维光子晶体能带的程序采用平面波传输法.docx

1、二维光子晶体能带的程序采用平面波传输法pwm法计算二维光子晶体能带的另一个程序已经分享了一个,这是另一个,编写方法不同,精度更高。程序完整含有注释,能同时显示TE和TM波的能带结构。略作修改也可分开显示。 已经分享了一个,这是另一个,编写方法不同,精度更高。程序完整含有注释,能同时显示TE和TM波的能带结构。略作修改也可分开显示。隐藏 function output_args = Untitled2( input_args )%UNTITLED2 Summary of this function goes here% Detailed explanation goes here%二维光子晶体带

2、结构计算%计算二维正方格子或三角格子，%介质圆柱（介电常数a）立于背景（介电常数b）之中clear; clc; tic; epssys=1.0e-6; %设定一个最小量，避免系统截断误差或除0错误%定义实际的正空间格子基矢latticetype=0; %定义格子类型，1 计算 正方格子，0 计算 三角格子switch latticetypecase 1a = 1;a1 = a* 1 0 0;a2 = a* 0 1 0;a3 = a* 0 0 1;case 0a = 1;a1 = a* sqrt(3)/2 1/2 0;a2 = a* -sqrt(3)/2 1/2 0;a3 = a* 0 0 1;

3、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);r

4、a2 = ra22(1:2);%选定参与运算的倒空间格矢量，即参与运算的平面波数量。%设定一个l,m的取值范围，变化l,m即可得出参与运算的平面波集合。NrSquare = 10; %选定k空间的尺度,即l,m（倒格矢G=l*b1+m*b2)的取值范围, NrSquare确定后，使用平面波数目可能为(2*NrSquare+1)2G = zeros(2*NrSquare+1)2,2); %初始化可能使用的倒格矢矩阵i = 1;for l = -NrSquare:NrSquarefor m = -NrSquare:NrSquareG(i,:) = l*ra1+m*ra2;i = i + 1;end

5、;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;% for l = -NrSquare:NrSquare% for m = -NrS

6、quare: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);for i=1:NGfor j=1:NGGij=norm(G(i,:)-G(j,:);if (Gij epssys)f(i,j)=(1/epsa)*Pf+(1/epsb)*(1-Pf);elsef(i,j)=(1/epsa-1/epsb)*Pf*2*bes

7、selj(1,Gij*Rc)/(Gij*Rc);end;end;end;%定义简约布里渊区的各高对称点switch latticetypecase 1T=(2*pi/a)*epssys 0;M=(2*pi/a)*1/2 1/2;X=(2*pi/a)*1/2 0;case 0T=(2*pi/a)*epssys 0;M=(2*pi/a)*sqrt(3)/3 0;K=(2*pi/a)*sqrt(3)/3 1/3;end;%对于简约布里渊区边界上的每个k，求解其特征频率THETA_TM=zeros(NG,NG); %待解的TM波矩阵THETA_TE=zeros(NG,NG); %待解的TE波矩阵Nkp

8、oints=10; %每个方向上取的点数，stepsize=0:1/(Nkpoints-1):1; %每个方向上的步长switch latticetypecase 1TX_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

9、 = zeros(Nkpoints,NG); %沿MT方向的TM波的待解的特征频率矩阵MT_TE_eig = zeros(Nkpoints,NG); %沿MT方向的TE波的待解的特征频率矩阵case 0TM_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波的

10、待解的特征频率矩阵KT_TM_eig = zeros(Nkpoints,NG); %沿KT方向的TM波的待解的特征频率矩阵KT_TE_eig = zeros(Nkpoints,NG); %沿KT方向的TE波的待解的特征频率矩阵end;for n=1:Nkpointsfprintf(n k-point:,int2str(n),of,int2str(Nkpoints),.n);%对于TX（正方格子）或TM（三角格子）方向上的每个k值，%求解其特征频率switch latticetypecase 1TX_step = stepsize(n)*(X-T)+T;case 0TM_step = steps

11、ize(n)*(M-T)+T;end;%先求非对角线上的元素for i=1:(NG-1)for j=(i+1):NGswitch latticetypecase 1kGi = TX_step+G(i,:);kGj = TX_step+G(j,:);case 0kGi = 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)=co

12、nj(THETA_TE(i,j);end;end;%再求对角线上的元素%for i=1:NGswitch latticetypecase 1kGi = TX_step+G(i,:);case 0kGi = 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矩阵的特征频率switch latticetypecase 1TX_TM_eig(n,:)=sort(sqrt(eig(THETA_T

13、M).;TX_TE_eig(n,:)=sort(sqrt(eig(THETA_TE).;case 0TM_TM_eig(n,:)=sort(sqrt(eig(THETA_TM).;TM_TE_eig(n,:)=sort(sqrt(eig(THETA_TE).;end;%对于XM（正方格子）或MK（三角格子）方向上的每个k值，%求解其特征频率switch latticetypecase 1XM_step = stepsize(n)*(M-X)+X;case 0MK_step = stepsize(n)*(K-M)+M;end;%先求非对角线上的元素for i=1:(NG-1)for j=(i+1

14、):NGswitch latticetypecase 1kGi = XM_step+G(i,:);kGj = XM_step+G(j,:);case 0kGi = 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;%再求对角线上的元素for i=1:NGswitch la

15、tticetypecase 1kGi = XM_step+G(i,:);case 0kGi = 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矩阵的特征频率switch latticetypecase 1XM_TM_eig(n,:)=sort(sqrt(eig(THETA_TM).;XM_TE_eig(n,:)=sort(sqrt(eig(THETA_TE).;case 0MK_TM

16、_eig(n,:)=sort(sqrt(eig(THETA_TM).;MK_TE_eig(n,:)=sort(sqrt(eig(THETA_TE).;end;%对于MT（正方格子）或KT（三角格子）方向上的每个k值，%求解其特征频率switch latticetypecase 1MT_step = stepsize(n)*(T-M)+M;case 0KT_step = stepsize(n)*(T-K)+K;end;%先求非对角线上的元素for i=1:(NG-1)for j=(i+1):NGswitch latticetypecase 1kGi = MT_step+G(i,:);kGj =

17、MT_step+G(j,:);case 0kGi = 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;%再求对角线上的元素for i=1:NGswitch latticetypecase 1kGi = MT_step+G(i,:);case 0kGi = KT_step

18、+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矩阵的特征频率switch latticetypecase 1MT_TM_eig(n,:)=sort(sqrt(eig(THETA_TM).;MT_TE_eig(n,:)=sort(sqrt(eig(THETA_TE).;case 0KT_TM_eig(n,:)=sort(sqrt(eig(THETA_TM).;KT_TE_eig(n,:)=sort(

19、sqrt(eig(THETA_TE).;end;end %这个是 for n=1:Nkpoints（115行）的循环结束fprintf(n Calculation Time: %d sec, toc)save pbs2D%以下部分是绘制光子晶体光子带图%首先将特定方向（正方格子:TX,XM,MT；三角格子：TM,MK,KT）离散化switch latticetypecase 1kaxis = 0;TXaxis = kaxis:norm(T-X)/(Nkpoints-1):(kaxis+norm(T-X);kaxis = kaxis + norm(T-X);XMaxis = kaxis:norm

20、(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);case 0kaxis = 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

21、 = 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)hold on;Nk=Nkpoints;for k=1:NGfor i=1:Nkpointsswitch latticetypecase 1EigFreq_TM(

22、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);case 0EigFreq_TM(i+0*Nk) = TM_TM_eig(i,k)

23、/(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; %for k=1:Nkpoints的结束switch latticetypecase 1plot(TXaxis(

24、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);case 0plot(TMaxis(1:Nk),EigFreq_TM(1+0*Nk:1*Nk),b,.MKaxis(1:N

25、k),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; %for k=1:NG的结束grid on;hold off;titlestr = 2D Photonic band structure for ;switch latticetypecas

26、e 1titlestr = titlestr, square ;case 0titlestr = titlestr, triangle ;end;titlestr = titlestr,lattice of dielectric columns.;titlestr = titlestr,(epsa=,num2str(epsa), epsb=, .num2str(epsb),);title (titlestr);xlabel(k-Space);ylabel(Frequency (omegaa/2piC);switch latticetypecase 1axis(0 MTaxis(Nkpoints) 0 1);set(gca,XTick,TXaxis(1).TXaxis(Nkpoints).XMaxis(Nkpoints).MTaxis(Nkpoints);xtixlabel = strvcat(T,X,M,T);case 0axis(0 KTaxis(Nkpoints) 0 1);set(gca,XTick,TMaxis(1).TMaxis(Nkpoints).MKaxis(Nkpoints).KTaxis(Nkpoints);xtixlabel = strvcat(T,M,K,T);end;set(gca,XTickLabel,xtixlabel);end