本程序为初学者使用,只考虑MT方向
下面的程序为matlab代码
只考虑MT方向
%This is a simple demo for Photonic Crystals simulation %This demo is for TE wave only, so only h wave is considered. %for TM direction only,10 points is considered. %---------------------------------------M %| / | %| / | %| / | %| --------------------|X %| T | %| | %| | %--------------------------------------- %equation :sum_{G‘,k}(K+G)(K+G‘)f(G-G‘)hz(k+G‘)=(omega/c)^2*hz(k+G) %G‘ can considerd as the index of column, and G as index of rows %[(K+G1)(K+G1)f(G1-G1) (K+G1)(K+G2)f(G1-G2) ][hz(G1)]=(omega/c)^2[hz(G1)] %[(K+G2)(K+G1)f(G2-G1) (K+G2)(K+G2)f(G2-G2) ][hz(G2)] [hz(G2)] %or: THETA_TE*Hz=(omega/c)^2*Hz %by Gao Haikuo %date:20170411 clear; clc; epssys=1.0e-6; %设定一个最小量,避免系统截断误差或除0错误 %this is the lattice vector and the reciprocal lattice vector a=1; a1=a*[1 0]; a2=a*[0 1]; b1=2*pi/a*[1 0];b2=2*pi/a*[0 1]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %定义晶格的参数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% epsa = 1; %介质柱的介电常数 epsb = 13; %背景的介电常数 Pf = 0.7; %Pf = Ac/Au 填充率,可根据需要自行设定 Au =a^2; %二维格子原胞面积 Rc = (Pf *Au/pi)^(1/2); %介质柱截面半径 Ac = pi*(Rc)^2; %介质柱横截面积 %construct the G list NrSquare = 10; NG =(2*NrSquare+1)^2; % NG is the number of the G value G = zeros(NG,2); i = 1; for l = -NrSquare:NrSquare for m = -NrSquare:NrSquare G(i,:)=l*b1+m*b2; i = i+1; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %生成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_TE=zeros(NG,NG); %待解的TE波矩阵 Nkpoints=10; %每个方向上取的点数, stepsize=0:1/(Nkpoints-1):1; %每个方向上的步长 TX_TE_eig = zeros(Nkpoints,NG); %沿MT方向的TE波的待解的特征频率矩阵 XM_TE_eig = zeros(Nkpoints,NG); %沿MT方向的TE波的待解的特征频率矩阵 MT_TE_eig = zeros(Nkpoints,NG); %沿MT方向的TE波的待解的特征频率矩阵 for n=1:Nkpoints %scan the 10 points along the TM direction fprintf([‘\n k-point:‘,int2str(n),‘of‘,int2str(Nkpoints),‘.\n‘]); MT_step = stepsize(n)*(T-M)+M; % get the k %先求非对角线上的元素 for i=1:(NG-1) % G for j=(i+1):NG % G‘ kGi = TX_step+G(i,:); %k+G kGj = TX_step+G(j,:); %K+G‘ THETA_TE(i,j)=f(i,j)*dot(kGi,kGj); %(K+G)(K+G‘)f(G-G‘) THETA_TE(j,i)=conj(THETA_TE(i,j)); end end %再求对角线上的元素 for i=1:NG kGi = TX_step+G(i,:); THETA_TE(i,i)=f(i,i)*norm(kGi)*norm(kGi); end MT_TE_eig(n,:)=sort(sqrt(eig(THETA_TE))).‘; end %draw 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; figure (1) hold on; Nk=Nkpoints; for k=1:NG for i=1:Nkpoints 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_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
时间: 2024-12-27 16:19:36