(PWEM前置知识)解析公式求相对介电常数分布函数的 Fourier 展开系数以及基于此的绘制实空间原胞的相对介电常数分布

正式计算前，检查下用傅里叶系数恢复的相对介电常数分布是很有必要的

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

================================================================
前面的现象可能找到了根本原因，就是反变换的公式估计理解错了。