11111.m

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

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
clear;
clear global;
close all;
eb=1;
ea=11.4;
ebxn=[ea,ea,ea,eb,ea;
      ea,ea,ea,eb,eb;
      eb,eb,ea,ea,ea;
      ea,eb,eb,ea,ea;
      ea,ea,eb,ea,ea];           
%The first quadrant of permittivity
pw=1;
cw=10;
Mx=7
My=7
flag=1

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);

        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
       
        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
     
        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);
disp('test ends');