收藏
下载资源下载资源 加入VIP,免费下载

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

上传人：b****3 文档编号：5510167 上传时间：2023-05-08 格式：DOCX 页数：15 大小：17.35KB
下载 相关 举报
二维光子晶体能带的程序采用平面波传输法.docx_第1页
第1页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第2页
第2页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第3页
第3页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第4页
第4页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第5页
第5页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第6页
第6页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第7页
第7页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第8页
第8页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第9页
第9页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第10页
第10页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第11页
第11页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第12页
第12页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第13页
第13页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第14页
第14页 / 共15页
二维光子晶体能带的程序采用平面波传输法.docx_第15页
第15页 / 共15页
亲，该文档总共15页，全部预览完了，如果喜欢就下载吧！
下载资源
资源描述

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

《二维光子晶体能带的程序采用平面波传输法.docx》由会员分享，可在线阅读，更多相关《二维光子晶体能带的程序采用平面波传输法.docx（15页珍藏版）》请在冰点文库上搜索。

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

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

pwm法计算二维光子晶体能带的另一个程序

已经分享了一个,这是另一个,编写方法不同,精度更高。

程序完整含有注释,能同时显示TE和TM波的能带结构。

略作修改也可分开显示。

已经分享了一个,这是另一个,编写方法不同,精度更高。

程序完整含有注释,能同时显示TE和TM波的能带结构。

略作修改也可分开显示。

<<隐藏

function[output_args]=Untitled2（input_args）

%UNTITLED2Summaryofthisfunctiongoeshere

%Detailedexplanationgoeshere

%二维光子晶体带结构计算

%计算二维正方格子或三角格子，

%介质圆柱（介电常数a）立于背景（介电常数b）之中

clear;clc;tic;epssys=1.0e-6;%设定一个最小量，避免系统截断误差或除0错误

%定义实际的正空间格子基矢

latticetype=0;%定义格子类型，1计算正方格子，0计算三角格子

switchlatticetype

case1

a=1;

a1=a*[100];

a2=a*[010];

a3=a*[001];

case0

a=1;

a1=a*[sqrt（3）/21/20];

a2=a*[-sqrt（3）/21/20];

a3=a*[001];

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;

forl=-NrSquare:

NrSquare

form=-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;

%forl=-NrSquare:

NrSquare

%form=-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）;

fori=1:

NG

forj=1:

NG

Gij=norm（G（i,:

）-G（j,:

））;

if（Gij

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;

%定义简约布里渊区的各高对称点

switchlatticetype

case1

T=（2*pi/a）*[epssys0];

M=（2*pi/a）*[1/21/2];

X=（2*pi/a）*[1/20];

case0

T=（2*pi/a）*[epssys0];

M=（2*pi/a）*[sqrt（3）/30];

K=（2*pi/a）*[sqrt（3）/31/3];

end;

%对于简约布里渊区边界上的每个k，求解其特征频率

THETA_TM=zeros（NG,NG）;%待解的TM波矩阵

THETA_TE=zeros（NG,NG）;%待解的TE波矩阵

Nkpoints=10;%每个方向上取的点数，

stepsize=0:

1/（Nkpoints-1）:

1;%每个方向上的步长

switchlatticetype

case1

TX_TM_eig=zeros（Nkpoints,NG）;%沿TX方向的TM波的待解的特征频率矩阵

TX_TE_eig=zeros（Nkpoints,NG）;%沿TX方向的TE波的待解的特征频率矩阵

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波的待解的特征频率矩阵

case0

TM_TM_eig=zeros（Nkpoints,NG）;%沿TM方向的TM波的待解的特征频率矩阵

TM_TE_eig=zeros（Nkpoints,NG）;%沿TM方向的TE波的待解的特征频率矩阵

MK_TM_eig=zeros（Nkpoints,NG）;%沿MK方向的TM波的待解的特征频率矩阵

MK_TE_eig=zeros（Nkpoints,NG）;%沿MK方向的TE波的待解的特征频率矩阵

KT_TM_eig=zeros（Nkpoints,NG）;%沿KT方向的TM波的待解的特征频率矩阵

KT_TE_eig=zeros（Nkpoints,NG）;%沿KT方向的TE波的待解的特征频率矩阵

end;

forn=1:

Nkpoints

fprintf（['\nk-point:

',int2str（n）,'of',int2str（Nkpoints）,'.\n']）;

%对于TX（正方格子）或TM（三角格子）方向上的每个k值，

%求解其特征频率

switchlatticetype

case1

TX_step=stepsize（n）*（X-T）+T;

case0

TM_step=stepsize（n）*（M-T）+T;

end;

%先求非对角线上的元素

fori=1:

（NG-1）

forj=（i+1）:

NG

switchlatticetype

case1

kGi=TX_step+G（i,:

）;

kGj=TX_step+G（j,:

）;

case0

kGi=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））;

THETA_TE（i,j）=f（i,j）*dot（kGi,kGj）;

THETA_TE（j,i）=conj（THETA_TE（i,j））;

end;

end;

%%%%%

%再求对角线上的元素

%%%%%

fori=1:

NG

switchlatticetype

case1

kGi=TX_step+G（i,:

）;

case0

kGi=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矩阵的特征频率

switchlatticetype

case1

TX_TM_eig（n,:

）=sort（sqrt（eig（THETA_TM）））.';

TX_TE_eig（n,:

）=sort（sqrt（eig（THETA_TE）））.';

case0

TM_TM_eig（n,:

）=sort（sqrt（eig（THETA_TM）））.';

TM_TE_eig（n,:

）=sort（sqrt（eig（THETA_TE）））.';

end;

%对于XM（正方格子）或MK（三角格子）方向上的每个k值，

%求解其特征频率

switchlatticetype

case1

XM_step=stepsize（n）*（M-X）+X;

case0

MK_step=stepsize（n）*（K-M）+M;

end;

%先求非对角线上的元素

fori=1:

（NG-1）

forj=（i+1）:

NG

switchlatticetype

case1

kGi=XM_step+G（i,:

）;

kGj=XM_step+G（j,:

）;

case0

kGi=MK_step+G（i,:

）;

kGj=MK_step+G（j,:

）;

end;

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;

%再求对角线上的元素

fori=1:

NG

switchlatticetype

case1

kGi=XM_step+G（i,:

）;

case0

kGi=MK_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;

%求解XM（正方格子）或MK（三角格子）方向上的k矩阵的特征频率

switchlatticetype

case1

XM_TM_eig（n,:

）=sort（sqrt（eig（THETA_TM）））.';

XM_TE_eig（n,:

）=sort（sqrt（eig（THETA_TE）））.';

case0

MK_TM_eig（n,:

）=sort（sqrt（eig（THETA_TM）））.';

MK_TE_eig（n,:

）=sort（sqrt（eig（THETA_TE）））.';

end;

%对于MT（正方格子）或KT（三角格子）方向上的每个k值，

%求解其特征频率

switchlatticetype

case1

MT_step=stepsize（n）*（T-M）+M;

case0

KT_step=stepsize（n）*（T-K）+K;

end;

%先求非对角线上的元素

fori=1:

（NG-1）

forj=（i+1）:

NG

switchlatticetype

case1

kGi=MT_step+G（i,:

）;

kGj=MT_step+G（j,:

）;

case0

kGi=KT_step+G（i,:

）;

kGj=KT_step+G（j,:

）;

end;

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;

%再求对角线上的元素

fori=1:

NG

switchlatticetype

case1

kGi=MT_step+G（i,:

）;

case0

kGi=KT_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矩阵的特征频率

switchlatticetype

case1

MT_TM_eig（n,:

）=sort（sqrt（eig（THETA_TM）））.';

MT_TE_eig（n,:

）=sort（sqrt（eig（THETA_TE）））.';

case0

KT_TM_eig（n,:

）=sort（sqrt（eig（THETA_TM）））.';

KT_TE_eig（n,:

）=sort（sqrt（eig（THETA_TE）））.';

end;

end%这个是'forn=1:

Nkpoints'（115行）的循环结束

fprintf（'\nCalculationTime:

%dsec',toc）

savepbs2D

%以下部分是绘制光子晶体光子带图

%首先将特定方向（正方格子:

TX,XM,MT；三角格子：

TM,MK,KT）离散化

switchlatticetype

case1

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）;

case0

kaxis=0;

TMaxis=kaxis:

norm（T-M）/（Nkpoints-1）:

（kaxis+norm（T-M））;

kaxis=kaxis+norm（T-M）;

MKaxis=kaxis:

norm（M-K）/（Nkpoints-1）:

（kaxis+norm（M-K））;

kaxis=kaxis+norm（M-K）;

KTaxis=kaxis:

norm（K-T）/（Nkpoints-1）:

（kaxis+norm（K-T））;

kaxis=kaxis+norm（K-T）;

end;

Ntraject=3;%所需绘制的特定方向的数目

EigFreq_TM=zeros（Ntraject*Nkpoints,1）;

EigFreq_TE=zeros（Ntraject*Nkpoints,1）;

figure

（1）

holdon;

Nk=Nkpoints;

fork=1:

NG

fori=1:

Nkpoints

switchlatticetype

case1

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）;

case0

EigFreq_TM（i+0*Nk）=TM_TM_eig（i,k）/（2*pi/a）;

EigFreq_TM（i+1*Nk）=MK_TM_eig（i,k）/（2*pi/a）;

EigFreq_TM（i+2*Nk）=KT_TM_eig（i,k）/（2*pi/a）;

EigFreq_TE（i+0*Nk）=TM_TE_eig（i,k）/（2*pi/a）;

EigFreq_TE（i+1*Nk）=MK_TE_eig（i,k）/（2*pi/a）;

EigFreq_TE（i+2*Nk）=KT_TE_eig（i,k）/（2*pi/a）;

end;

end;%fork=1:

Nkpoints的结束

switchlatticetype

case1

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'）;

case0

plot（TMaxis（1:

Nk）,EigFreq_TM（1+0*Nk:

1*Nk）,'b',...

MKaxis（1:

Nk）,EigFreq_TM（1+1*Nk:

2*Nk）,'b',...

KTaxis（1:

Nk）,EigFreq_TM（1+2*Nk:

3*Nk）,'b'）;

plot（TMaxis（1:

Nk）,EigFreq_TE（1+0*Nk:

1*Nk）,'r',...

MKaxis（1:

Nk）,EigFreq_TE（1+1*Nk:

2*Nk）,'r',...

KTaxis（1:

Nk）,EigFreq_TE（1+2*Nk:

3*Nk）,'r'）;

end;

end;%fork=1:

NG的结束

gridon;

holdoff;

titlestr='2DPhotonicbandstructurefor';

switchlatticetype

case1

titlestr=[titlestr,'square'];

case0

titlestr=[titlestr,'triangle'];

end;

titlestr=[titlestr,'latticeofdielectriccolumns.'];

titlestr=[titlestr,'（epsa=',num2str（epsa）,'epsb=',...

num2str（epsb）,'）'];

title（titlestr）;

xlabel（'k-Space'）;

ylabel（'Frequency（\omegaa/2\piC）'）;

switchlatticetype

case1

axis（[0MTaxis（Nkpoints）01]）;

set（gca,'XTick',[TXaxis

（1）...

TXaxis（Nkpoints）...

XMaxis（Nkpoints）...

MTaxis（Nkpoints）]）;

xtixlabel=strvcat（'T','X','M','T'）;

case0

axis（[0KTaxis（Nkpoints）01]）;

set（gca,'XTick',[TMaxis

（1）...

TMaxis（Nkpoints）...

MKaxis（Nkpoints）...

KTaxis（Nkpoints）]）;

xtixlabel=strvcat（'T','M','K','T'）;

end;

set（gca,'XTickLabel',xtixlabel）;

end

展开阅读全文
相关资源
猜你喜欢
相关搜索
资源标签

当前位置：首页 > IT计算机 > 电脑基础知识

copyright@ 2008-2023 冰点文库 网站版权所有

经营许可证编号:鄂ICP备19020893号-2