squarepixel_ebxn_2d_fc.asv

通过计算晶体结构来确定其光谱分布

function ret=SquarePixel_ebxn_2D_FC(ebxn,pw,Mx,My,flag)
%return Fourier Coefficients of 2D ebxn with Square pixels
%ebxn=the first quadrant in permittivity 
%pw=pixel width
%cw=cell width
%Mx,My=Half width 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
chN=t(1);
cw=2*chN*pw;        %Number of pixels
Nx=2*Mx+1;
Ny=2*My+1;
N=Nx*Ny;
dl=(-2*Mx:2*Mx);
dm=(-2*My:2*My);
ret=zeros(N);
switch flag
    case 1
        tx=zeros(chN,Nx,Nx);
        t=zeros(1,4*Mx+1);
        ebxn1=1./ebxn;
        for y=1:chN
            for l=1:4*Mx+1
                t(l)=SquarePixel_ebxn_1D_FC(ebxn1(y,:),pw,dl(l));
            end
            tx(y,:,:)=inv(toeplitz(t(Nx:4*Mx+1),fliplr(t(1:Nx))));
        end
        clear t;
        t=zeros(4*My+1,Nx,Nx);
        for l=1:Nx
            for l1=1:Nx
                for m=1:4*My+1
                    t(m,l,l1)=SquarePixel_ebxn_1D_FC((tx(:,l,l1)).',pw,dm(m));
                end
            end
        end
        clear tx;
        tr=zeros(N,Nx);
        tc=zeros(N,Nx);
        for m=1:Ny
            tc((m-1)*Nx+1:m*Nx,:)=t(Ny-1+m,:,:);
            tr((m-1)*Nx+1:m*Nx,:)=t(Ny+1-m,:,:);
        end
        ret=bktoeplitz(tc,tr);
        clear t tr tc;
    case 2
        ty=zeros(chN,Ny,Ny);
        t=zeros(1,4*My+1);        
        ebxn1=1./ebxn;
        for x=1:chN
            for l=1:4*My+1
                t(l)=SquarePixel_ebxn_1D_FC((ebxn1(:,x)).',pw,dm(l));
            end
            ty(x,:,:)=inv(toeplitz(t(Ny:4*My+1),fliplr(t(1:Ny))));
        end
        clear t;
        t=zeros(4*Mx+1,Ny,Ny);
        for m=1:Ny
            for m1=1:Ny
                for l=1:4*Mx+1
                    t(l,m,m1)=SquarePixel_ebxn_1D_FC((ty(:,m,m1)).',pw,dl(l));
                end
            end
        end
        clear ty;
        tr=zeros(N,Ny);
        tc=zeros(N,Ny);
        for l=1:Nx
            tc((l-1)*Ny+1:l*Ny,:)=t(Nx-1+l,:,:);
            tr((l-1)*Ny+1:l*Ny,:)=t(Nx+1-l,:,:);
        end
        ret=bktoeplitz(tc,tr);
        clear t tr tc;
    case 3
        tx=zeros(chN,Nx,Nx);
        t=zeros(1,4*Mx+1);
        ebxn1=ebxn;
        for y=1:chN
            for l=1:4*Mx+1
                t(l)=SquarePixel_ebxn_1D_FC(ebxn1(y,:),pw,dl(l));
            end
            tx(y,:,:)=toeplitz(t(Nx:4*Mx+1),fliplr(t(1:Nx)));
        end
        clear t;
        t=zeros(4*My+1,Nx,Nx);
        for l=1:Nx
            for l1=1:Nx
                for m=1:4*My+1
                    t(m,l,l1)=SquarePixel_ebxn_1D_FC((tx(:,l,l1)).',pw,dm(m));
                end
            end
        end
        clear tx;
        tr=zeros(N,Nx);
        tc=zeros(N,Nx);
        for m=1:Ny
            tr((m-1)*Nx+1:m*Nx,:)=t(Ny-1+m,:,:);
            tc((m-1)*Nx+1:m*Nx,:)=t(Ny+1-m,:,:);
        end
        ret=bktoeplitz(tc,tr);
        clear t tr tc;
end