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

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; %设定一个最小量，避免系统截断误差或除0错误 %定义实际的正空间格子基矢 latticetype=0; %定义格子类型，1 计算 正方格子，0 计算 三角格子 switch latticetype case 1 a = 1; a1 = a* [1 0 0]; a2 = a* [0 1 0]; a3 = a* [0 0 1]; case 0 a = 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/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; for l = -NrSquare:NrSquare for m = -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; % for l = -NrSquare:NrSquare % for m = -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); 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; %定义简约布里渊区的各高对称点 switch latticetype case 1 T=(2*pi/a)*[epssys 0]; M=(2*pi/a)*[1/2 1/2]; X=(2*pi/a)*[1/2 0]; case 0 T=(2*pi/a)*[epssys 0];