• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号

二维光子晶体带隙仿真Matlab完全程序_平面波展开法

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

本程序为二维光子晶体Matlab仿真程序,该结果与文献【1】Molding the flow of light,p68 figure 2相互吻合

 

 

主程序

%This is a simple demo for Photonic Crystals simulation
%10 points is considered.
%---------------------------------------M
%| / |
%| / |
%| / |
%| --------------------|X
%| T |
%| |
%| |
%---------------------------------------
%by Gao Haikuo
%date:20170411
clear; clc; epssys=1.0e-6; %设定一个最小量,避免系统截断误差或除0错误
global NG G f Nkpoints
%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];
Nkpoints=10; %每个方向上取的点数,
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%定义晶格的参数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
epsa = 8.9; %inner
epsb = 1; %outer
Pf = 0.1257; %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];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[MT_TE_eig]=getEigValue(M,T,0);
[TX_TE_eig]=getEigValue(T,X,0);
[XM_TE_eig]=getEigValue(X,M,0);
[MT_TM_eig]=getEigValue(M,T,1);
[TX_TM_eig]=getEigValue(T,X,1);
[XM_TM_eig]=getEigValue(X,M,1);
%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);
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);
end
linehandle=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',...
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',...
'LineWidth',1 );
set (linehandle, 'linesmoothing', 'on');
end
grid on;
xlabel('K-Space');
yLabel('Frequency(\omegaa/2\piC)');
axis([0 MTaxis(Nkpoints) 0 0.8]);
set(gca,'XTick',[TXaxis(1), TXaxis(Nkpoints),...
XMaxis(Nkpoints),MTaxis(Nkpoints)]);
xtixlabel = strvcat('T','X','M','T');
set(gca,'XTickLabel',xtixlabel);

核心获取本征值的函数

function [eigValue]=getEigValue(k_begin,k_end,mode)
%mode :0 for TE 1 for TH
global NG G f Nkpoints
THETA=zeros(NG,NG); %待解的TE波矩阵
stepsize=0:1/(Nkpoints-1):1; %每个方向上的步长
TX_TE_eig = zeros(Nkpoints,NG); 
for n=1:Nkpoints %scan the 10 points along the direction
    fprintf(['\n k-point:',int2str(n),'of',int2str(Nkpoints),'.\n']);
    step = stepsize(n)*(k_end-k_begin)+k_begin;  % get the k
    for i=1:(NG-1)   % G
        for j=(i+1):NG % G'
            kGi = step+G(i,:); %k+G
            kGj = step+G(j,:); %K+G'
            switch mode
                case 0  %TE mode
                    THETA(i,j)=f(i,j)*dot(kGi,kGj); %(K+G)(K+G')f(G-G')
                case 1 %TH mode
                     THETA(i,j)=f(i,j)*norm(kGi)*norm(kGj);
            end
            THETA(j,i)=conj(THETA(i,j)); 
        end
    end    
    for i=1:NG
        kGi = step+G(i,:);
        THETA(i,i)=f(i,i)*norm(kGi)*norm(kGi); 
    end
    eigValue(n,:)=sort(sqrt(eig(THETA))).';
end

 


鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
【转】Matlab混合编程之NET组件(C#篇)发布时间:2022-07-22
下一篇:
MATLAB之易经卜卦程序+GUI发布时间:2022-07-22
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap