ebxn_2d_fs.m

来自「通过计算晶体结构来确定其光谱分布」· M 代码 · 共 52 行

function ret=ebxn_2D_FS(ebxn,pw,cw,m,n,flag)
%return Fourier series of 2D ebxn
%ebxn=the first quadrant in permittivity 
%pw=pixel width
%cw=cell width
%m,n=orders index of Fourier coefficient along x and y axis,respectively
%flag=1: inverse rule for x axis, Laurent's rule for y axis;
%flag=2: inverse rule for y axis, Laurent's rule for x axis;
%flag=3: Laurent's rule for both x and y axes
t=size(ebxn);
if t(1)~=t(2)
    ret=-1;
    disp('ebxn is not a square matrix');
    return;
end
M=t(1);
t0=4*pw^2/cw^2;
mt=pi*m*pw/cw;
nt=pi*n*pw/cw;
if mt==0
    sm=1;
else
    sm=sin(mt)/mt;
end
nt=pi*n*pw/cw;
if nt==0
    sn=1;
else
    sn=sin(nt)/nt;
end
et=zeros(1,M);
switch flag
    case 1 
        ei=1./ebxn;
        for y=1:M
            et(y)=1/((sum(ei(:,y).*cos(mt*(2*(1:M)'-1)))*sm)+eps);
        end
        ret=t0*sum(et.*cos(nt*(2*(1:M)-1)))*sn;
    case 2
        ei=1./ebxn;
        for x=1:M
            et(x)=1/((sum(ei(x,:).*cos(nt*(2*(1:M)-1)))*sn)+eps);
        end
        ret=t0*sum(et.*cos(mt*(2*(1:M)-1)))*sm;
    case 3
        for y=1:M
            et(y)=sum(ebxn(:,y).*cos(mt*(2*(1:M)'-1)))*sm;
        end
        ret=t0*sum(et.*cos(nt*(2*(1:M)-1)))*sn;
end