二维体光子晶体的平面波展开法代码
2024-06-23 06:47:18
%书上的代码，和FEM符合的更好
%在这个代码里试着把单位原胞的相对介电常数分布画出来
%这个代码的单位原胞的中心就是(0,0)点，也就是坐标原点
%The program for the computation of the PhC photonic
%band structure for 2D PhC.
%Parameters of the structure are defined by the PhC
%period, elements radius, and by the permittivities
%of elements and background material.
%Input parameters: PhC period, radius of an
%element, permittivities
%Output data: The band structure of the 2D PhC.
clear all;clc;close all
%The variable a defines the period of the structure.
%It influences on results in case of the frequency
%normalization.
a=1e-6; % in meters，晶格常数，单位m
%The variable r contains elements radius. Here it is
%defined as a part of the period.
r=0.4*a; % in meters，单位m，介质柱半径
%The variable eps1 contains the information about
%the relative background material permittivity.
eps1=9;%背景材料介电常数
%The variable eps2 contains the information about
%permittivity of the elements composing PhC. In our
%case permittivities are set to form the structure
%perforated membrane
eps2=1;%介质柱的相对介电常数
%The variable precis defines the number of k-vector
%points between high symmetry points
precis=20;%和BZ区的离散有关
%The variable nG defines the number of plane waves.
%Since the number of plane waves cannot be
%arbitrary value and should be zero-symmetric, the
%total number of plane waves may be determined
%as (nG*2-1)^2
nG=15;%和平面波的数量有关
%The variable precisStruct defines contains the
%number of nodes of discretization mesh inside in
%unit cell discretization mesh elements.
precisStruct=50;%单位原胞离散份数，最好设为偶数
%The following loop carries out the definition
%of the unit cell. The definition is being
%made in discreet form by setting values of inversed
%dielectric function to mesh nodes.
nx=1;
for countX=-a/2:a/precisStruct:a/2%就是网格节点的x坐标
ny=1;
for countY=-a/2:a/precisStruct:a/2%就是网格节点的y坐标
%The following condition allows to define the circle
%with of radius r
if(sqrt(countX^2+countY^2)<r)%判断网格节点坐标是否处于中心的介质柱内
%Setting the value of the inversed dielectric
%function to the mesh node
struct(nx,ny)=1/eps2;
%Saving the node coordinate
xSet(nx)=countX;
ySet(ny)=countY;
else
struct(nx,ny)=1/eps1;
xSet(nx)=countX;
ySet(ny)=countY;
end
ny=ny+1;
end
nx=nx+1;
end
%The computation of the area of the mesh cell. It is
%necessary for further computation of the Fourier
%expansion coefficients.
dS=(a/precisStruct)^2;%单位微元的面积
%Forming 2D arrays of nods coordinates
xMesh=meshgrid(xSet(1:length(xSet)-1));
yMesh=meshgrid(ySet(1:length(ySet)-1))';
%Transforming values of inversed dielectric function
%for the convenience of further computation
structMesh=struct(1:length(xSet)-1,...
1:length(ySet)-1)*dS/(max(xSet)-min(xSet))^2;%就是(C10)的一部分%Defining the k-path within the Brilluoin zone. The
%k-path includes all high symmetry points and precis
%points between them
kx(1:precis+1)=0:pi/a/precis:pi/a;
ky(1:precis+1)=zeros(1,precis+1);
kx(precis+2:precis+precis+1)=pi/a;
ky(precis+2:precis+precis+1)=...
pi/a/precis:pi/a/precis:pi/a;
kx(precis+2+precis:precis+precis+1+precis)=...
pi/a-pi/a/precis:-pi/a/precis:0;
ky(precis+2+precis:precis+precis+1+precis)=...
pi/a-pi/a/precis:-pi/a/precis:0;
%After the following loop, the variable numG will
%contain real number of plane waves used in the
%Fourier expansion
numG=1;
%The following loop forms the set of reciprocal
%lattice vectors.
for Gx=-nG*2*pi/a:2*pi/a:nG*2*pi/a
for Gy=-nG*2*pi/a:2*pi/a:nG*2*pi/a
G(numG,1)=Gx;
G(numG,2)=Gy;
numG=numG+1;
end
end
%numG为82，实际上多了1
%The next loop computes the Fourier expansion
%coefficients which will be used for matrix
%differential operator computation.
for countG=1:numG-1Mark0=countG/(numG-1)
for countG1=1:numG-1temp=exp(1i*((G(countG,1)-G(countG1,1))*...
xMesh+(G(countG,2)-G(countG1,2))*yMesh));CN2D_N(countG,countG1)=sum(sum(structMesh.*...
temp));
end
end%The next loop computes matrix differential operator
%in case of TE mode. The computation
%is carried out for each of earlier defined
%wave vectors.
for countG=1:numG-1Mark=countG/(numG-1)
for countG1=1:numG-1
for countK=1:length(kx)
M(countK,countG,countG1)=...
CN2D_N(countG,countG1)*((kx(countK)+G(countG,1))*...
(kx(countK)+G(countG1,1))+...
(ky(countK)+G(countG,2))*(ky(countK)+G(countG1,2)));
end
end
end
%The computation of eigen-states is also carried
%out for all wave vectors in the k-path.
for countK=1:length(kx)
%Taking the matrix differential operator for current
%wave vector.
MM(:,:)=M(countK,:,:);
%Computing the eigen-vectors and eigen-states of
%the matrix
[D,V]=eig(MM);
%Transforming matrix eigen-states to the form of
%normalized frequency.
dispe(:,countK)=sqrt(V*ones(length(V),1))*a/2/pi;
end
%Plotting the band structure
%Creating the output field
figure(3);
%Creating axes for the band structure output.
ax1=axes;
%Setting the option
%drawing without cleaning the plot
hold on;
%Plotting the first 8 bands
for u=1:8
plot(abs(dispe(u,:)),'r','LineWidth',2);
%If there are the PBG, mark it with blue rectangle
if(min(dispe(u+1,:))>max(dispe(u,:)))
rectangle('Position',[1,max(dispe(u,:)),...
length(kx)-1,min(dispe(u+1,:))-...
max(dispe(u,:))],'FaceColor','b',...
'EdgeColor','b');
end
end
%Signing Labeling the axes
%下面是和FEM算的数据比较，没有的话可以删除下面这段%%
Data=csvread('C:\Users\LUO1\Desktop\ComsolDate.csv');
uniX=unique(Data(:,1));
Ymat=[];
for zz=1:length(uniX)Ymat=[Ymat,Data(Data(:,1)==uniX(zz),2)];
end
plot(1:size(Ymat,2),Ymat,'bo','MarkerSize',3)
%%%%%绘制FEM的数据%%%%%%%set(ax1,'xtick',...
[1 precis+1 2*precis+1 3*precis+1]);
set(ax1,'xticklabel',['G';'X';'M';'G']);
ylabel('Frequency \omegaa/2\pic','FontSize',16);
xlabel('Wavevector','FontSize',16);
xlim([1 length(kx)])
set(ax1,'XGrid','on');title('book')
蓝色小圆是FEM绘制的能带。我原先按照自己的理解改动了下上面的代码，但是和FEM的结果符合得没上面的代码好。
抄的书上的代码: Photonic Crystals Physics and Practical Modeling (Igor A. Sukhoivanov, Igor V. Guryev (auth.))
下面是相同参数，但按照我自己的理解改了一点点的代码，结果和FEM的结果符合得不是很好，有无这方面的大佬评论区解释下原因？或者任何见解也行
%这份代码和FEM符合的不是很好
%在这个代码里试着把单位原胞的相对介电常数分布画出来
%这个代码的单位原胞的中心就是(0,0)点，也就是坐标原点
%The program for the computation of the PhC photonic
%band structure for 2D PhC.
%Parameters of the structure are defined by the PhC
%period, elements radius, and by the permittivities
%of elements and background material.
%Input parameters: PhC period, radius of an
%element, permittivities
%Output data: The band structure of the 2D PhC.
clear all;clc;close all
%The variable a defines the period of the structure.
%It influences on results in case of the frequency
%normalization.
a=1e-6; % in meters，晶格常数，单位m
%The variable r contains elements radius. Here it is
%defined as a part of the period.
r=0.4*a; % in meters，单位m，介质柱半径
%The variable eps1 contains the information about
%the relative background material permittivity.
eps1=9;%背景材料介电常数
%The variable eps2 contains the information about
%permittivity of the elements composing PhC. In our
%case permittivities are set to form the structure
%perforated membrane
eps2=1;%介质柱的相对介电常数
%The variable precis defines the number of k-vector
%points between high symmetry points
precis=20;%和BZ区的离散有关
%The variable nG defines the number of plane waves.
%Since the number of plane waves cannot be
%arbitrary value and should be zero-symmetric, the
%total number of plane waves may be determined
%as (nG*2-1)^2
nG=15;%和平面波的数量有关
%The variable precisStruct defines contains the
%number of nodes of discretization mesh inside in
%unit cell discretization mesh elements.
precisStruct=50;%单位原胞离散份数，最好设为偶数
%The following loop carries out the definition
%of the unit cell. The definition is being
%made in discreet form by setting values of inversed
%dielectric function to mesh nodes.
nx=1;
for countX=-a/2:a/precisStruct:a/2%就是网格节点的x坐标
ny=1;
for countY=-a/2:a/precisStruct:a/2%就是网格节点的y坐标
%The following condition allows to define the circle
%with of radius r
if(sqrt(countX^2+countY^2)<r)%判断网格节点坐标是否处于中心的介质柱内
%Setting the value of the inversed dielectric
%function to the mesh node
struct(nx,ny)=1/eps2;
%Saving the node coordinate
xSet(nx)=countX;
ySet(ny)=countY;
else
struct(nx,ny)=1/eps1;
xSet(nx)=countX;
ySet(ny)=countY;
end
ny=ny+1;
end
nx=nx+1;
end
%The computation of the area of the mesh cell. It is
%necessary for further computation of the Fourier
%expansion coefficients.
dS=(a/precisStruct)^2;%单位微元的面积
%Forming 2D arrays of nods coordinates
xMesh=meshgrid(xSet(1:length(xSet)));
yMesh=meshgrid(ySet(1:length(ySet)))';
%Transforming values of inversed dielectric function
%for the convenience of further computation
structMesh=struct(1:length(xSet),...
1:length(ySet))*dS/(max(xSet)-min(xSet))^2;%就是(C10)的一部分%Defining the k-path within the Brilluoin zone. The
%k-path includes all high symmetry points and precis
%points between them
kx(1:precis+1)=0:pi/a/precis:pi/a;
ky(1:precis+1)=zeros(1,precis+1);
kx(precis+2:precis+precis+1)=pi/a;
ky(precis+2:precis+precis+1)=...
pi/a/precis:pi/a/precis:pi/a;
kx(precis+2+precis:precis+precis+1+precis)=...
pi/a-pi/a/precis:-pi/a/precis:0;
ky(precis+2+precis:precis+precis+1+precis)=...
pi/a-pi/a/precis:-pi/a/precis:0;
%After the following loop, the variable numG will
%contain real number of plane waves used in the
%Fourier expansion
numG=1;
%The following loop forms the set of reciprocal
%lattice vectors.
for Gx=-nG*2*pi/a:2*pi/a:nG*2*pi/a
for Gy=-nG*2*pi/a:2*pi/a:nG*2*pi/a
G(numG,1)=Gx;
G(numG,2)=Gy;
numG=numG+1;
end
end
%numG为82，实际上多了1
%The next loop computes the Fourier expansion
%coefficients which will be used for matrix
%differential operator computation.
for countG=1:numG-1Mark0=countG/(numG-1)
for countG1=1:numG-1temp=exp(1i*((G(countG,1)-G(countG1,1))*...
xMesh+(G(countG,2)-G(countG1,2))*yMesh));CN2D_N(countG,countG1)=sum(sum(structMesh.*...
temp));
end
end%The next loop computes matrix differential operator
%in case of TE mode. The computation
%is carried out for each of earlier defined
%wave vectors.
for countG=1:numG-1Mark=countG/(numG-1)
for countG1=1:numG-1
for countK=1:length(kx)
M(countK,countG,countG1)=...
CN2D_N(countG,countG1)*((kx(countK)+G(countG,1))*...
(kx(countK)+G(countG1,1))+...
(ky(countK)+G(countG,2))*(ky(countK)+G(countG1,2)));
end
end
end
%The computation of eigen-states is also carried
%out for all wave vectors in the k-path.
for countK=1:length(kx)
%Taking the matrix differential operator for current
%wave vector.
MM(:,:)=M(countK,:,:);
%Computing the eigen-vectors and eigen-states of
%the matrix
[D,V]=eig(MM);
%Transforming matrix eigen-states to the form of
%normalized frequency.
dispe(:,countK)=sqrt(V*ones(length(V),1))*a/2/pi;
end
%Plotting the band structure
%Creating the output field
figure(3);
%Creating axes for the band structure output.
ax1=axes;
%Setting the option
%drawing without cleaning the plot
hold on;
%Plotting the first 8 bands
for u=1:8
plot(abs(dispe(u,:)),'r','LineWidth',2);
%If there are the PBG, mark it with blue rectangle
if(min(dispe(u+1,:))>max(dispe(u,:)))
rectangle('Position',[1,max(dispe(u,:)),...
length(kx)-1,min(dispe(u+1,:))-...
max(dispe(u,:))],'FaceColor','b',...
'EdgeColor','b');
end
end
Data=csvread('C:\Users\LUO1\Desktop\ComsolDate.csv');
uniX=unique(Data(:,1));
Ymat=[];
for zz=1:length(uniX)Ymat=[Ymat,Data(Data(:,1)==uniX(zz),2)];
end
plot(1:size(Ymat,2),Ymat,'bo','MarkerSize',3)
%Signing Labeling the axes
set(ax1,'xtick',...
[1 precis+1 2*precis+1 3*precis+1]);
set(ax1,'xticklabel',['G';'X';'M';'G']);
ylabel('Frequency \omegaa/2\pic','FontSize',16);
xlabel('Wavevector','FontSize',16);
xlim([1 length(kx)])
set(ax1,'XGrid','on');title('test')和FEM的对比：
%================================================%
我明白为什么我改动后的代码是错的了，原因是需要把xMesh、yMesh当作每个离散小微元的中心坐标，所以书上的代码更准确。
上面的理解又是错的，书上代码参考的那个数值求傅里叶系数的公式应该理解成二维离散傅里叶级数(2D DFS)，总之这个知识点的却非常细微，不注意很容易以为自己搞懂了，但实际上理解是错的。正确的离散方法求介电分布函数的傅里叶系数的公式如下：
这实际上是个2D DFS，但fft2也可以间接计算这个式子。
上图是正确的正空间原胞划分网格示意图，四个红点勾勒出了正空间原胞的大小，蓝色的点就是计算中用到的节点数据。中间的红点是(0,0)点，红圈是介质柱轮廓。
=======================================
%=======================================%
为何下面两段代码的结果之差不是严格等于0？非常奇怪。有无这方面的大佬解释下？
Kappa(countG,countG1)*...(norm([kx(countK)+G(countG,1),ky(countK)+G(countG,2)])*...norm([kx(countK)+G(countG1,1),ky(countK)+G(countG1,2)]))
Kappa(countG,countG1)*...
norm([kx(countK)+G(countG,1),ky(countK)+G(countG,2)])*...
norm([kx(countK)+G(countG1,1),ky(countK)+G(countG1,2)])
写PWEM算法的时候只有用第一段代码得到的结果才是对的。
%=======================================%
%=========================================%
在某些nG参数下，书上代码绘制的能带会全部挤到0附近不知道为何，
但是如果加大precisStruct结果又会变回正常；在使用数值方法计算相对介电函数的傅里叶系数时nG和precisStruct参数的相对大小需要小心设置
后面我找到了这一现象可能的原因，和离散傅里叶变换(DFT)，离散傅里叶级数(DFS)的公式具体形式有关。基于此我猜测有关系：2nG+1应小于或等于precisStruct（本质是2nG+1大小应小于等于正空间原胞离散份数）
%=========================================%
对书上代码的改写(TE、TM都可以计算，数值或解析方法算相对介电函数的傅里叶系数)
实际上无论算TE 还是TM能带 数值方法算的相对介电函数的傅里叶系数得到的能带反而比解析方法得到的能带和FEM方法误差更小（实际上解析方法算的相对介电函数比数值方法的质量好些），不清楚为何会这样。
图像如下：
代码：
%为什么按照我的想法去设置微元的中心坐标效果不行？
%数值法算Kappa的结果和FEM更接近，不知道为何
clear ;clc;close all
TEorTM='TM';
NUMorANA='ANA';%数值或解析方法算Kappa矩阵(傅里叶系数)
a=1e-6; % in meters，晶格常数，单位mr=0.4*a; % in meters，单位m，介质柱半径eps1=9;%背景材料介电常数eps2=1;%介质柱的相对介电常数precis=20;%和BZ区的离散有关nG=16;%和平面波的数量有关precisStruct=44;%单位原胞离散段数(x、y方向)，最好设为偶数，
%因此节点数为precisStruct+1%%%正空间基矢(可自己设置)%%%
a1=[a,0,0];%为方便求倒格子基矢，必须是三维矢量
a2=[0,a,0];
%%%正空间基矢%%%%%%%%%%%%%求倒格子基矢%%%%%
a3=[0,0,1];%单位向量
Omega=abs(a3*cross(a1,a2)');%正空间单位原胞面积
s=2*pi./Omega;
b1(1,1)=a2(1,2);b1(1,2)=-a2(1,1);
b2(1,1)=-a1(1,2);b2(1,2)=a1(1,1);
b1=b1.*s;%行数组，最终结果(必须是行数组)
b2=b2.*s;%行数组，最终结果(必须是行数组)
%%%%求倒格子基矢%%%%%nodeSet=-a/2:a/precisStruct:a/2;
len_nodeSet=length(nodeSet);
%初始化
struct=zeros(len_nodeSet,len_nodeSet);
xMesh0=zeros(len_nodeSet,len_nodeSet);
yMesh0=zeros(len_nodeSet,len_nodeSet);for nx=1:len_nodeSetcountX=nodeSet(nx);for ny=1:len_nodeSetcountY=nodeSet(ny);if(sqrt(countX^2+countY^2)<r)%判断网格节点坐标是否处于中心的介质柱内struct(nx,ny)=1/eps2;%其实就是每个节点的相对介电常数xMesh0(nx,ny)=countX;%其实就是每个节点的x坐标yMesh0(nx,ny)=countY;%其实就是每个节点的y坐标elsestruct(nx,ny)=1/eps1;xMesh0(nx,ny)=countX;yMesh0(nx,ny)=countY;endendend
xMesh0(:,end)=[];
xMesh0(end,:)=[];
yMesh0(:,end)=[];
yMesh0(end,:)=[];%其实xSet=ySet=[-a/2:a/precisStruct:a/2]dS=(a/precisStruct)^2;%单位微元的面积
%Forming 2D arrays of nods coordinates%注意xMesh=xMesh0.'
%注意yMesh=yMesh0.'
%注意xMesh=yMesh0
%注意yMesh=xMesh0
%xMesh0、yMesh0其实是节点的坐标？
%=======================%
%=======================%
figure
for ii = 1:numel(xMesh0)plot(xMesh0(ii),yMesh0(ii),'b.','Markersize',8)hold onrectangle('Position',[xMesh0(ii)-0.5*(a/precisStruct),...yMesh0(ii)-0.5*(a/precisStruct),(a/precisStruct),(a/precisStruct)]);end
xx1=a/2;
siteSto=[xx1,xx1;xx1,-xx1;-xx1,xx1;-xx1,-xx1];
hold on
scatter(siteSto(:,1),siteSto(:,2),'r.')
hold on
plot(0,0,'r.','Markersize',8)
hold on
% x y是中心，r是半径
rectangle('Position',[0-r,0-r,2*r,2*r],'Curvature',[1,1],'edgecolor',[1,0,0])
hold off
axis equal
title('xMesh、yMesh')
%=======================%
%=======================%structMesh=struct(1:precisStruct,...1:precisStruct)*dS/Omega;%就是公式(C10)的一部分
%注意struct(1:length(xSet)-1,1:length(ySet)-1)等价于
%struct(:,end)=[];struct(end,:)=[];%========离散不可约布里渊区(只适用于正方晶格)==========%
kx(1:precis+1)=0:pi/a/precis:pi/a;
ky(1:precis+1)=zeros(1,precis+1);
kx(precis+2:precis+precis+1)=pi/a;
ky(precis+2:precis+precis+1)=...pi/a/precis:pi/a/precis:pi/a;
kx(precis+2+precis:precis+precis+1+precis)=...pi/a-pi/a/precis:-pi/a/precis:0;
ky(precis+2+precis:precis+precis+1+precis)=...pi/a-pi/a/precis:-pi/a/precis:0;%========离散不可约布里渊区(只适用于正方晶格)==========%%The following loop forms the set of reciprocal
%lattice vectors.
G=zeros((2*nG+1)^2,2);
numG=0;
for m=-nG:nGfor n=-nG:nGnumG=numG+1;G(numG,1:2)=m.*b1+n.*b2;end
end
%G里面没有重复的G_{||}向量
numG=size(G,1);
%The next loop computes the Fourier expansion
%coefficients which will be used for matrix
%differential operator computation.
%遍历G数组Kappa=zeros(numG,numG);
Indexset=-nG:nG;
if strcmp(NUMorANA,'NUM')==1%数值法计算Kappafor countG=1:numGMark0=countG/(numG)for countG1=1:numGtemp=exp(-1i*((G(countG,1)-G(countG1,1)).*...xMesh0+(G(countG,2)-G(countG1,2)).*yMesh0));%矩阵Kappa(countG,countG1)=sum(sum(structMesh.*...temp));%其实就是kappa(G_{||}-G_{||}^{\prime}})构成的矩阵%G(countG,1:2)其实就是G_{||}，%G(countG1,1:2)其实就是G_{||}^{\prime}endend%%%%%%%%计算kappa(G_{||})%%%%%%%%Gmat=zeros(length(Indexset),length(Indexset));Kapmat=zeros(length(Indexset),length(Indexset));for m=1:length(Indexset)for n=1:length(Indexset)g0=Indexset(m).*b1+Indexset(n).*b2;%行数组Gmat(m,n)=g0(1)+1j.*g0(2);temp=exp(-1i*(g0(1).*...xMesh0+g0(2).*yMesh0));%矩阵Kapmat(m,n)=sum(sum(structMesh.*...temp));endend%%%%%%%%%%%%%%%%%%%%%%%endf=(pi*r^2)/Omega;
if strcmp(NUMorANA,'ANA')==1%解析法计算Kappafor countG=1:numGMark0=countG/(numG)for countG1=1:numGtemp=norm(G(countG,1:2)-G(countG1,1:2))*r;%|G|rif temp==0Kappa(countG,countG1)=1/eps1+f*(1/eps2-1/eps1);elseKappa(countG,countG1)=2*f*(1/eps2-1/eps1)*besselj(1,temp)/temp;end%其实就是kappa(G_{||}-G_{||}^{\prime}})构成的矩阵%G(countG,1:2)其实就是G_{||}，%G(countG1,1:2)其实就是G_{||}^{\prime}endend%%%%%%%%计算kappa(G_{||})%%%%%%%%Gmat=zeros(length(Indexset),length(Indexset));Kapmat=zeros(length(Indexset),length(Indexset));for m=1:length(Indexset)for n=1:length(Indexset)g0=Indexset(m).*b1+Indexset(n).*b2;%行数组normG=norm(g0);%|G|Gmat(m,n)=g0(1)+1j.*g0(2);if Indexset(m)==0 && Indexset(n)==0Kapmat(m,n)=1/eps1+f*(1/eps2-1/eps1);elseKapmat(m,n)=2*f*(1/eps2-1/eps1)*besselj(1,(normG*r))/(normG*r);endendend%%%%%%%%%%%%%%%%%%%%%%%end%======绘制Kappa(G_{||})还原的相对介电常数分布==========%
%G_{||}-G_{||}^{\prime}}包含了很多重复元素
%因此kappa(G_{||}-G_{||}^{\prime}})不能用来恢复相对介电常数
%分布,只能用kappa(G_{||})
EpsMat=zeros(precisStruct+1,precisStruct+1);
for uu=1:precisStruct+1uu/(precisStruct+1)X=nodeSet(uu);for vv=1:precisStruct+1Y=nodeSet(vv);EpsMat(uu,vv)=sum(sum(Kapmat.*exp(1j*(X.*real(Gmat)+Y.*imag(Gmat)))));end
end
figure
mesh(1./real(EpsMat)) %画图
title([NUMorANA,' RealEps',' e1=',num2str(eps1),' e2=',num2str(eps2)])figure
mesh(1./imag(EpsMat)) %画图
title([NUMorANA,' ImagEps',' e1=',num2str(eps1),' e2=',num2str(eps2)])figure
imagesc(1./real(EpsMat))
colorbar%======绘制Kappa(G_{||})还原的相对介电常数分布==========%%构建特征矩阵M=zeros(numG,numG,length(kx));%特征矩阵初始化
for countK=1:length(kx)Mark=countK/length(kx)for countG=1:numGfor countG1=1:numGif strcmp(TEorTM,'TE')==1%TEM(countG,countG1,countK)=...Kappa(countG,countG1)*...([kx(countK)+G(countG,1),ky(countK)+G(countG,2)]*...([kx(countK)+G(countG1,1),ky(countK)+G(countG1,2)].'));endif strcmp(TEorTM,'TM')==1%TMM(countG,countG1,countK)=...Kappa(countG,countG1)*...(norm([kx(countK)+G(countG,1),ky(countK)+G(countG,2)])*...norm([kx(countK)+G(countG1,1),ky(countK)+G(countG1,2)]));endendend
end%The computation of eigen-states is also carried
%out for all wave vectors in the k-path.
for countK=1:length(kx)%Taking the matrix differential operator for current%wave vector.MM(:,:)=M(:,:,countK);%Computing the eigen-vectors and eigen-states of%the matrix[D,V]=eig(MM);%Transforming matrix eigen-states to the form of%normalized frequency.dispe(:,countK)=sqrt(V*ones(length(V),1))*a/2/pi;
end
%Plotting the band structure
%Creating the output field
figure;
%Creating axes for the band structure output.
ax1=axes;
%Setting the option
%drawing without cleaning the plot
hold on;
%Plotting the first 8 bands
for u=1:8plot(abs(dispe(u,:)),'r','LineWidth',2);%If there are the PBG, mark it with blue rectangleif(min(dispe(u+1,:))>max(dispe(u,:)))rectangle('Position',[1,max(dispe(u,:)),...length(kx)-1,min(dispe(u+1,:))-...max(dispe(u,:))],'FaceColor','b',...'EdgeColor','b');end
end
%Signing Labeling the axes
%下面是和FEM算的数据比较，没有的话可以删除下面这段%%
if strcmp(TEorTM,'TE')==1Data=csvread('C:\Users\LUO1\Desktop\ComsolDateTE.csv');
end
if strcmp(TEorTM,'TM')==1Data=csvread('C:\Users\LUO1\Desktop\ComsolDateTM.csv');
enduniX=unique(Data(:,1));
Ymat=[];
for zz=1:length(uniX)Ymat=[Ymat,Data(Data(:,1)==uniX(zz),2)];
endplot(1:size(Ymat,2),Ymat,'bo','MarkerSize',3)
hold off
%%%%%绘制FEM的数据%%%%%%%set(ax1,'xtick',...[1 precis+1 2*precis+1 3*precis+1]);
set(ax1,'xticklabel',['G';'X';'M';'G']);
ylabel('Frequency \omegaa/2\pic','FontSize',16);
xlabel('Wavevector','FontSize',16);
xlim([1 length(kx)])
set(ax1,'XGrid','on');
title([TEorTM,' ',NUMorANA,' book'])实际上nG不是越大越收敛于FEM的能带结果，它和precisStruct大小有个很有趣的关系。
precisStruct最好设为偶数(经验发现，设为偶数，TE模式的收敛效果会更好)，2nG+1<=precisStruct(按照DFS的定义)，并且2nG+1不应和precisStruct相差很大。
二维体光子晶体的平面波展开法代码相关推荐
- matlab在光子晶体,利用平面波展开法在matlab中计算一维光子晶体的带隙结构
黄晓泽 周琦 饶黄云 摘 要:利用平面波展开法,在王建东所研究的光在光子晶体(photonic crystal)中传输所满足的非线性薛定谔方程基础上,计算光子晶体间隙结构的本征方程,根据本征方程在m ...
- 光子晶体中的平面波展开法学习
把平面波展开法当做一个响应器,即有输入和输出,先谈它的输入,输入有: 1.介质柱的相对介电常数 2.背景的相对介电常数 3.填充率(介质柱横截面积比上正格子中二维原胞面积,三角格子光子晶体的原胞是一个 ...
- matlab晶体能带,matlab平面波展开法的二维光子晶体能带研究+程序
摘 要 :二维光子晶体可以作为对光子传输控制的新型材料.本文主要通过平面波展开法对二维光子晶体进行数值计算及其性质分析.首先我们介绍了二维光子晶体的基础概念.结构.介电性能等特性.然后基于麦克斯韦方 ...
- 二维光子晶体禁带的遗传优化算法实现
二维光子晶体禁带的遗传优化算法MATLAB源代码 光子晶体中因周期性结构而存在的频率禁带称为光子禁带,光子禁带的存在是光子晶体具有广泛应用前景的重要原因. 禁带越大,可控光的频带也越宽,因此如何设计合 ...
- OpenCV学习笔记（二十六）——小试SVM算法ml OpenCV学习笔记（二十七）——基于级联分类器的目标检测objdect OpenCV学习笔记（二十八）——光流法对运动目标跟踪Video Ope
OpenCV学习笔记(二十六)--小试SVM算法ml 总感觉自己停留在码农的初级阶段,要想更上一层,就得静下心来,好好研究一下算法的东西.OpenCV作为一个计算机视觉的开源库,肯定不会只停留在数字图 ...
- 【光学】基于Matlab实现二维光子晶体的能带图和场
1 内容介绍 为了了解电磁波在光子晶体中的传输特性,用MATLAB与时域有限差分法把电磁波在真空与光子晶体中的传播实时可视化,并给出了场的空间静态分布.数值模拟的结果表明,禁带中的波被光子晶体控制,其 ...
- 机器学习（二）：k近邻法（kNN）
引言 一.k近邻法的三要素 1.距离度量 2.k值的选择 3.分类决策规则 二.k近邻算法及代码实现(python) 1.算法 2.代码实现(python) 三.k近邻法的实现:kd树 1.构造kd树 ...
- 一阶指数平滑c语言,时间序列数据之一、二、三阶指数平滑法
原标题:时间序列数据之一.二.三阶指数平滑法 基本概念 时间序列数据: 是指一个数据序列,在时域上存在相互依赖的关系 前提假设: 时间序列分析一般假设我们获得的数据在时域上具有一定的相互依赖关系,例如 ...
- 【数据结构与算法】3.二叉搜索树(BST)、c++代码
二叉搜索树(BST).c++代码 参考 https://blog.csdn.net/c_living/article/details/81021510
最新文章
- BMC：幼年特发性关节炎患儿肠道菌群的特征、生物标记的识别及其在临床预测中的作用...
- python正确的赋值语句-Python 赋值语句技巧之序列解包
- linux httpd 域名映射,51CTO博客-专业IT技术博客创作平台-技术成就梦想
- spark rdd Transformation和Action 剖析
- 阿里技术专家推荐的20本书，免费送！
- Android之提示Unable to instantiate fragment***MyLikeFragment .could not find Fragment constructor
- mysql的indata文件_【数据分析】MySQL之不能导入本地文件“Loading local data is disable;”...
- Huaman Gene Functions
- java flash 压缩_Java和flash通信中数据的zlib压缩与解压缩
- 每天学一点flash(6) FLASH 8 和 FLASH CS3 加载外部文本的区别 (转载)
- 分片报文的最后一片_分片报文攻击防范 - CloudEngine 12800, 12800E V200R005C10 配置指南-安全 - 华为...
- 为什么我的 WordPress 网站被封了？
- 【直击DTCC】浪潮霍俊路详解数据库集群锁机制管理
- Tcp/ip协议 详解
- PHP之50个开源项目
- KVM远程迁移启动报错
- IntelliJ IDEA设置编码格式
- TexStudio 快捷键 ctrl+R 查找替换
- 手机里tencent文件夹能删吗_【我要删除tencent文件夹】手机删了tencent文件夹_tencent什么文件夹...
- NOIP模拟赛csy2021/10/30