piecewise_band_1d.m

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

function [OH,OE]=piecewise_band_1D(ebxn,pw,M,ON,Nk,kx,kz,Gz,k0)
%return Fourier Coefficients of 1D ebxn for H and E-polarization
%ebxn=the first quadrant in permittivity 
%pw=pixel width
%M=Half width of Fourier coefficient along z axis
%k0=modulus of maximal k
t=size(ebxn);
chN=max(t(1));
cw=2*chN*pw;        %Number of pixels
N=2*M+1;
dl=(-2*M:2*M);
clear t;
t=zeros(1,4*M+1);
ebxn1=1./ebxn;
for l=1:4*M+1
    t(l)=SquarePixel_ebxn_1D_FC(ebxn1,pw,dl(l));
end
FCz=inv(toeplitz(t(N:4*M+1),fliplr(t(1:N))));
t=zeros(1,4*M+1);
ebxn1=ebxn;
for l=1:4*M+1
    t(l)=SquarePixel_ebxn_1D_FC(ebxn1,pw,dl(l));
end
FCx=toeplitz(t(N:4*M+1),fliplr(t(1:N)));

omegaH=zeros(Nk,N);
omegaE=zeros(Nk,N);

h=waitbar(0,'Eigen value...');
for m=1:Nk
    waitbar(m/Nk,h);
    Kx=eye(N)*kx(m)+Gz;
    Kz=eye(N)*kz(m)+Gz;
    MH=Kx*inv(FCz)*Kx+Kz*inv(FCx)*Kz;
    ME=inv(FCx)*(Kx.^2+Kz.^2);
    [V,D]=eig(MH);
    for n=1:N
        omegaH(m,n)=D(n,n);
    end
    [V,D]=eig(ME);    
    for n=1:N
        omegaE(m,n)=D(n,n);
    end
end
close(h);
for m=1:Nk
    omegaH(m,:)=sort(omegaH(m,:));
    omegaE(m,:)=sort(omegaE(m,:));
end
OHt=zeros(Nk,ON);
OEt=zeros(Nk,ON);
for n=1:Nk
    for m=1:N
        if omegaH(n,m)>=0
            mt=N-m+1;
            if mt>ON
                OHt(n,:)=omegaH(n,m:m+ON-1);
            else
                OHt(n,1:mt)=omegaH(n,m:m+mt-1);
            end
            break;
        end
    end
end
for n=1:Nk
    for m=1:N
        if omegaE(n,m)>=0
            mt=N-m+1;
            if mt>ON
                OEt(n,:)=omegaE(n,m:m+ON-1);
            else
                OEt(n,1:mt)=omegaE(n,m:m+mt-1);
            end
            break;
        end
    end
end
OH=sqrt(OHt)/k0;
OE=sqrt(OEt)/k0;