raw-resourse.zip_二维光子晶体_二维声子晶体_光子晶体带隙_声子晶体带隙_带隙图

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %二维光子晶体带结构计算 %计算二维正方格子或三角格子， %介质圆柱（介电常数a）立于背景（介电常数b）之中 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; clc; tic; epssys=1.0e-6; %设定一个最小量，避免系统截断误差或除0错误 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %定义实际的正空间格子基矢 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% a = 1; a1 = a* [1 0 0]; a2 = a* [0 1 0]; a3 = a* [0 0 1]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %定义晶格的参数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% epsa = 1; %介质柱的介电常数 epsb = 13; %背景的介电常数 Pf = 0.7; %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)); 19 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,:); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %生成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; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %定义简约布里渊区的各高对称点 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% T=(2*pi/a)*[epssys 0]; M=(2*pi/a)*[1/2 1/2]; X=(2*pi/a)*[1/2 0]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对于简约布里渊区边界上的每个k，求解其特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% THETA_TM=zeros(NG,NG); %待解的TM波矩阵 THETA_TE=zeros(NG,NG); %待解的TE波矩阵 Nkpoints=10; %每个方向上取的点数， stepsize=0:1/(Nkpoints-1):1; %每个方向上的步长 TX_TM_eig = zeros(Nkpoints,NG); %沿TX方向的TM波的待解的特征频率矩阵 TX_TE_eig = zeros(Nkpoints,NG); %沿TX方向的TE波的待解的特征频率矩阵 21 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波的待解的特征频率矩阵 for n=1:Nkpoints fprintf(['\n k-point:',int2str(n),'of',int2str(Nkpoints),'.\n']); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对于TX（正方格子）或TM（三角格子）方向上的每个k值， %求解其特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TX_step = stepsize(n)*(X-T)+T; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %先求非对角线上的元素 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:(NG-1) for j=(i+1):NG kGi = TX_step+G(i,:); kGj = TX_step+G(j,:); 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:NG kGi = TX_step+G(i,:); 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矩阵的特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TX_TM_eig(n,:)=sort(sqrt(eig(THETA_TM))).'; TX_TE_eig(n,:)=sort(sqrt(eig(THETA_TE))).'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对于XM（正方格子）或MK（三角格子）方向上的每个k值， %求解其特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% XM_step = stepsize(n)*(M-X)+X; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %先求非对角线上的元素 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:(NG-1) for j=(i+1):NG kGi = XM_step+G(i,:); kGj = XM_step+G(j,:); 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:NG kGi = XM_step+G(i,:); 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矩阵的特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% XM_TM_eig(n,:)=sort(sqrt(eig(THETA_TM))).'; XM_TE_eig(n,:)=sort(sqrt(eig(THETA_TE))).'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对于MT（正方格子）或KT（三角格子）方向上的每个k值， %求解其特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MT_step = stepsize(n)*(T-M)+M; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %先求非对角线上的元素 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:(NG-1) for j=(i+1):NG kGi = MT_step+G(i,:); kGj = MT_step+G(j,:); 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:NG kGi = MT_step+G(i,:); THETA_TM(i,i)=f(i,i)*norm(kGi)*norm(kGi); THETA_TE(i,i)=f(i,i)*norm(kGi)*norm(kGi); end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 25 %求解TX（正方格子）或TM（三角格子）方向上的k矩阵的特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MT_TM_eig(n,:)=sort(sqrt(eig(THETA_TM))).'; MT_TE_eig(n,:)=sort(sqrt(eig(THETA_TE))).'; end fprintf('\n Calculation Time: %d sec', toc) save pbs2D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %以下部分是绘制光子晶体光子带图 %首先将特定方向（正方格子:TX,XM,MT；三角格子：TM,MK,KT）离散化 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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); Ntraject = 3; %所需绘制的特定方向的数目 EigFreq_TM=zeros(Ntraject*Nkpoints,1); EigFreq_TE=zeros(Ntraject*Nkpoints,1); figure (1) hold on; Nk=Nkpoints; for k=1:NG for i=1:Nkpoints 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); end; 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'); end; grid on; hold off; titlestr = '2D Photonic band structure for '; titlestr = [titlestr, ' square ']; titlestr = [titlestr,'lattice of dielectric columns.']; titlestr = [titlestr,'(epsa=',num2str(epsa), ' epsb=', ... num2str(epsb),')']; title (titlestr); xlabel('k-Space'); ylabel('Frequency (\omegaa/2\piC)'); axis([0 MTaxis(Nkpoints) 0 1]); set(gca,'XTick',[TXaxis(1)... TXaxis(Nkpoints)... XMaxis(Nkpoints)... MTaxis(Nkpoints)]); xtixlabel = strvcat('T','X','M','T'); set(gca,'XTickLabel',xtixlabel);