正式计算前,检查下用傅里叶系数恢复的相对介电常数分布是很有必要的
function []=Mian_Fourier()%解析公式求相对介电常数分布函数的 Fourier 展开系数(也就是笔记中的\kappa(G_{||})),这一步是很快的
%注意是正方晶格,正空间单位原胞是矩形的!!!坐标原点在单位原胞几何中心
clear;clc;close all
a=1;%晶格常数(复式格子的情况可能要斟酌一下)
xlength=a;%单位原胞沿x方向的长度
ylength=a;%单位原胞沿y方向的长度
Omega=xlength*ylength;%正空间单位原胞面积
numNode=61;%原胞离散网格的节点数(最好设置为奇数整数个)
Order=60;%倒格矢G的展开阶数
f=0.5;%填充比,就是pi*(ra^2)/Omega,大小不要离0或1太近
ra=sqrt(f*Omega/pi);%介质柱的半径
ea=3;%无限高介质柱的相对介电常数
eb=1;%背景材料的相对介电常数
%背景材料和柱子的相对介电常数不能相差太大,要不然会导致傅里叶系数失效
tem_eps=(1/ea)-(1/eb);%一个参数;
if mod(numNode,2)==0%numNode偶数
disp('警告: 网格离散节点numNode设置为了偶数,可能会出现错误!')
end
%%%%求倒格子基矢%%%%%
a1=[xlength,0,0];%正空间基矢
a2=[0,ylength,0];%三维
a3=[0,0,1];%单位向量
A=abs(a3*cross(a1,a2)');
s=2*pi./A;
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;%行数组,最终结果(必须是行数组)
%%%%求倒格子基矢%%%%%
kapMat=zeros(2*Order+1,2*Order+1);%kappa(G_{||})矩阵初始化
mnSto=[];
for m=-Order:Order
for n=-Order:Order
G=m.*b1+n.*b2;%倒格子矢G_{||},行数组
normG=sqrt(G(1)^2+G(2)^2);%|G_{||}|
if m==0&&n==0
kapMat(m+Order+1,n+Order+1)=1/eb+f*tem_eps;
else
kapMat(m+Order+1,n+Order+1)=2*f*tem_eps*besselj(1,normG*ra)/(normG*ra);
end
end
end
%计算结束,开始绘图
EpsDraw(kapMat,numNode,xlength,ylength,b1,b2,Order)
end
function []=EpsDraw(kapMat,numNode,xlength,ylength,b1,b2,Order)%绘制实空间原胞的相对介电常数分布,这一步很慢
EpsMat=zeros(numNode,numNode);
numNode2=floor(numNode/2);%floor(X) 将X的每个元素舍入为小于或等于该元素的最近整数
rxSto=zeros(numNode,1);%列数组
rySto=zeros(numNode,1);%列数组
for uu=-numNode2:numNode2
rxSto(uu+numNode2+1,1)=uu*xlength/(numNode-1);%存储了正空间网格节点的实际x坐标,列数组
end
for vv=-numNode2:numNode2
rySto(vv+numNode2+1,1)=vv*ylength/(numNode-1);%存储了正空间网格节点的实际y坐标,列数组
end
count=0;
for uu=-numNode2:numNode2%和实空间坐标有关
count=count+1;
Mark=count/(numNode)
for vv=-numNode2:numNode2%和实空间坐标有关
rx=rxSto(uu+numNode2+1,1);
ry=rySto(vv+numNode2+1,1);
for m=-Order:Order%和G_{||}有关
for n=-Order:Order%和G_{||}有关
G=m.*b1+n.*b2;%倒格子矢G_{||},行数组
Gx=G(1);
Gy=G(2);
EpsMat(uu+numNode2+1,vv+numNode2+1)=EpsMat(uu+numNode2+1,vv+numNode2+1)+...
kapMat(m+Order+1,n+Order+1)*exp(1j*(Gx*rx+Gy*ry));
end
end
end
end
figure(1)
mesh(real(1./EpsMat)) %画图
title('RealEps')
figure(2)
mesh(imag(1./EpsMat)) %画图
title('ImagEps')
figure(3)
imagesc(real(1./EpsMat))
colorbar
end
上述代码运用的公式:
这都是些常见公式。
注意里面的Mian_Fourier函数计算是极快的,慢的是EpsDraw部分,在真正用PWEM计算能带之前很有必要检查下Mian_Fourier计算得到的傅里叶展开系数在实空间的相对介电常数分布图像(也就是EpsDraw的工作),因为很有可能在设定的参数下该分布图像完全不准确,就像下面这种情况一样。
里面的ea、eb参数不能相差过大,比如ea=13,eb=1,就会出现原胞中相对介电常数分布失真的情况,如下图所示(边缘部分误差极大):
这种情况下需要将代码中的Order参数设得很大才会有较好的收敛效果。
这个情况目前还不知道如何解决(23.5.3),留到以后吧。
====================================================
找到了一个这种震荡的可能原因—傅里叶级数分解中的吉布斯现象(https://www.zhihu.com/question/68135290),这一现象是所有基于傅里叶级数的方法的通病。
=====================================================
使用了数值方法代替解析方法算傅里叶系数还是会有尖锐的峰,并且相同参数下数值算法的峰更尖锐
=====================================================
上述代码的一个快速版本
function []=Mian_Fourier()%解析公式求相对介电常数分布函数的 Fourier 展开系数(也就是笔记中的\kappa(G_{||})),这一步是很快的
%注意是正方晶格,正空间单位原胞是矩形的!!!坐标原点在单位原胞几何中心
clear;clc;close all
a=1;%晶格常数(复式格子的情况可能要斟酌一下)
xlength=a;%单位原胞沿x方向的长度
ylength=a;%单位原胞沿y方向的长度
Omega=xlength*ylength;%正空间单位原胞面积
numNode=61;%原胞离散网格的节点数(最好设置为奇数整数个)
Order=60;%倒格矢G的展开阶数
f=0.5;%填充比,就是pi*(ra^2)/Omega,大小不要离0或1太近
ra=sqrt(f*Omega/pi);%介质柱的半径
ea=1;%无限高介质柱的相对介电常数
eb=9;%背景材料的相对介电常数
%背景材料和柱子的相对介电常数不能相差太大,要不然会导致傅里叶系数失效
tem_eps=(1/ea)-(1/eb);%一个参数;
if mod(numNode,2)==0%numNode偶数
disp('警告: 网格离散节点numNode设置为了偶数,可能会出现错误!')
end
%%%%求倒格子基矢%%%%%
a1=[xlength,0,0];%正空间基矢
a2=[0,ylength,0];%三维
a3=[0,0,1];%单位向量
A=abs(a3*cross(a1,a2)');
s=2*pi./A;
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;%行数组,最终结果(必须是行数组)
%%%%求倒格子基矢%%%%%
kapMat=zeros(2*Order+1,2*Order+1);%kappa(G_{||})矩阵初始化
Gmat=kapMat;
mnSto=[];
for m=-Order:Order
for n=-Order:Order
G=m.*b1+n.*b2;%倒格子矢G_{||},行数组
normG=sqrt(G(1)^2+G(2)^2);%|G_{||}|
if m==0&&n==0
kapMat(m+Order+1,n+Order+1)=1/eb+f*tem_eps;
else
kapMat(m+Order+1,n+Order+1)=2*f*tem_eps*besselj(1,normG*ra)/(normG*ra);
end
Gmat(m+Order+1,n+Order+1)=G(1)+1j*G(2);
end
end
%计算结束,开始绘图
EpsDraw(kapMat,numNode,xlength,ylength,b1,b2,Order,Gmat)
end
function []=EpsDraw(kapMat,numNode,xlength,ylength,b1,b2,Order,Gmat)%绘制实空间原胞的相对介电常数分布,这一步很慢
EpsMat=zeros(numNode,numNode);
numNode2=floor(numNode/2);%floor(X) 将X的每个元素舍入为小于或等于该元素的最近整数
rxSto=zeros(numNode,1);%列数组
rySto=zeros(numNode,1);%列数组
for uu=-numNode2:numNode2
rxSto(uu+numNode2+1,1)=uu*xlength/(numNode-1);%存储了正空间网格节点的实际x坐标,列数组
end
for vv=-numNode2:numNode2
rySto(vv+numNode2+1,1)=vv*ylength/(numNode-1);%存储了正空间网格节点的实际y坐标,列数组
end
count=0;
for uu=-numNode2:numNode2%和实空间坐标有关
count=count+1;
Mark=count/(numNode)
for vv=-numNode2:numNode2%和实空间坐标有关
rx=rxSto(uu+numNode2+1,1);
ry=rySto(vv+numNode2+1,1);
Gx=real(Gmat);%矩阵
Gy=imag(Gmat);
EpsMat(uu+numNode2+1,vv+numNode2+1)=sum(sum(kapMat.*exp(1j*(rx.*Gx+ry.*Gy))));
end
end
figure(1)
mesh(real(1./EpsMat)) %画图
title('RealEps')
figure(2)
mesh(imag(1./EpsMat)) %画图
title('ImagEps')
figure(3)
imagesc(real(1./EpsMat))
colorbar
end
================================================================
前面的现象可能找到了根本原因,就是反变换的公式估计理解错了。