二维光子晶体能带的程序采用平面波传输法精品资料

1、pwm 法计算二维光子晶体能带的另一个程序已经分享了一个,这是另一个 , 编写方法不同 ,精度更高。程序完整含有注释,能同时显示TE 和 TM 波的能带结构。略作修改也可分开显示。function output_args = Untitled2( input_args )%UNTITLED2 Summary of this function goes here% Detailed explanation goes here% 二维光子晶体带结构计算%计算二维正方格子或三角格子，%介质圆柱（介电常数a）立于背景（介电常数b）之中clear; clc; tic; epssys=1.0e-6; %

2、设定一个最小量，避免系统截断误差或除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;end;%定义晶格的参数epsa = 1; %介质柱的介电常数epsb = 10.5; %背景的介电常数Pf = 0.28274; %Pf = Ac/A

3、u填充率，可根据需要自行设定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 即可得出参与运算的平面波

4、集合。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;end;NG = i - 1; % 实际使用的平面波数目G = G(1:NG ,:);%还有另一种选取的原则，即考虑到园

5、对称的话，则尽可能使平面波波矢的集合均匀分布%在一个球内，根据这样的原则，可先选取一个比较大的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 = -NrSquare:NrSquare% Glm = l*ra1+m*ra2;% if (norm(Glm<=Gm

6、ax)% 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*besselj(1,Gij*Rc)/(Gij*Rc);end;end;end;%定义

7、简约布里渊区的各高对称点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 波矩阵Nkpoints=10; % 每个方向上取的点数，stepsize=

8、0:1/(Nkpoints-1):1; % 每个方向上的步长switch latticetypecase 1TX_TM_eig = zeros(Nkpoints,NG); %TX_TE_eig = zeros(Nkpoints,NG); %XM_TM_eig = zeros(Nkpoints,NG); %XM_TE_eig = zeros(Nkpoints,NG); %MT_TM_eig = zeros(Nkpoints,NG); %MT_TE_eig = zeros(Nkpoints,NG); %沿 TX 方向的 TM 波的待解的特征频率矩阵沿 TX 方向的 TE 波的待解的特征频率矩阵沿

9、XM 方向的 TM 波的待解的特征频率矩阵沿 XM 方向的 TE 波的待解的特征频率矩阵沿 MT 方向的 TM 波的待解的特征频率矩阵沿 MT 方向的 TE 波的待解的特征频率矩阵case 0TM_TM_eig = zeros(Nkpoints,NG); % TM_TE_eig = zeros(Nkpoints,NG); %沿 沿TM TM方向的方向的TM 波的待解的特征频率矩阵TE 波的待解的特征频率矩阵MK_TM_eig = zeros(Nkpoints,NG); %沿 MK 方向的 TM 波的待解的特征频率矩阵MK_TE_eig = zeros(Nkpoints,NG); %沿 MK 方

10、向的 TE 波的待解的特征频率矩阵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 1T

11、X_step = stepsize(n)*(X-T)+T;case 0TM_step = stepsize(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);

12、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 = 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

13、 latticetypecase 1TX_TM_eig(n,:)=sort(sqrt(eig(THETA_TM).'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):NGswitch latticetypecase 1kGi = XM_step+G(i,:);kGj =