squarepixelband_2d.m

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

%SquarePixelBand_2D.m

ebxn0=8.8542e-12;
mu0=4*pi*1e-7;
c=1/sqrt(ebxn0*mu0);
t=size(ebxn);
cw=t(1)*pw*2;
%Cell width
t0=cat(2,fliplr(ebxn),ebxn);
ebxn0=cat(1,t0,flipud(t0));

hN=0;
hN=hN+1;
figure(hN);
eM=max(max(ebxn0));
em=min(min(ebxn0));
sM=64;
sm=0;
r=(sM-sm)/(eM-em);
%set(gca,'YLim',[0 10])
%set(gca,'XLim',[0 10])
ttt=(ebxn0-em)*r+sm;
xx=find(ttt==0);
ttt(xx)=63;
ttt(find(ttt==64))=0;

image(ttt);
%colormap(gray);
cc=[0,0,1;1,1,1];
colormap(cc);
grid on
title('Original data');

%Maximal half width of spectrum
Nx=2*Mx+1;
Ny=2*My+1;
N=Nx*Ny;
%Step number for calculation of k  

FCz=SquarePixel_ebxn_2D_FC(ebxn,pw,Mx,My,3);
if FCz==-1
    return;
end
FCx=SquarePixel_ebxn_2D_FC(ebxn,pw,Mx,My,1);
if FCx==-1
    return;
end
FCy=SquarePixel_ebxn_2D_FC(ebxn,pw,Mx,My,2);
if FCy==-1
    return;
end
gx0=(-Mx:Mx)*2*pi/cw;
gy0=(-My:My)*2*pi/cw;
gx=zeros(1,N);
gy=zeros(1,N);
for n=1:2*My+1
    for m=1:2*Mx+1
        gx(m+(n-1)*Nx)=gx0(m);
        gy(m+(n-1)*Ny)=gy0(n);
    end
end
Gx=diag(gx,0);
Gy=diag(gy,0);

Nk=3*Nk0;
k0=2*pi/cw;
k1=pi/cw;
kx0=0;
kx1=k1;
ky0=0;
ky1=0;
kx=linspace(0,k1,Nk0);
ky=linspace(0,0,Nk0);
kx=cat(2,kx,linspace(k1,k1,Nk0));
ky=cat(2,ky,linspace(0,k1,Nk0));
kx=cat(2,kx,linspace(k1,0,Nk0));
ky=cat(2,ky,linspace(k1,0,Nk0));
FCx=real(FCx);
xx=find(abs(FCx)<10e-13);
FCx(xx)=0;
FCy=real(FCy);
yy=find(abs(FCy)<10e-13);
FCy(yy)=0;
FCz=real(FCz);
zz=find(abs(FCz)<10e-13);
FCz(zz)=0;
disp('start compute eigen')

str={'Eigen value...'};
[OH,OE]=piecewise_band_2D(FCx,FCy,FCz,Gx,Gy,ON,kx,ky,k0,Nk,N,str);

findbandgap;
pictureofgap;

disp('SquarePixelBand ends');