这个程序结合超胞能带计算技巧，对局域态密度和电子密度进行可视化

nx=30;  % 
ny=8; % 体系宽度（y方向的长度）
[x,y]=zigzag_graphene(nx,ny);
N=4*nx*ny; % 原子总数
xx=x(N/3+1:N*2/3,1);  % 取出中间device区域的坐标
yy=y(N/3+1:N*2/3,1);  % 取出中间device区域的坐标
%plot(x,y,'.','MarkerSize',20)
t1=-2.7;
t2=0.0038/3/sqrt(3);
H=Hamiltonian_NN_graphene(x,y,t1);
Hsp=Hamiltonian_Haldane(x,y,sqrt(3),t2);
H=H+Hsp;
N=length(H);
HDL=H(N/3+1:N*2/3,1:N/3);
HD=H(N/3+1:N*2/3,N/3+1:N*2/3);
HDR=H(N/3+1:N*2/3,N*2/3+1:N);
n = length(HD);
kx=3;
Hk=HDL*exp(-1i*kx*nx/3)+HD+HDR*exp(1i*kx*nx/3);
[V,E]=eig(Hk);
plot(diag(E))
ans=V(:,n/2);
local_density = V(:,n/2).*conj(V(:,n/2));
m = length(xx);
Zs=0;  % 总的局域电子密度
sigma=0.5;  % 高斯函数的标准差
du=0.01; % x方向采样间隔
dv=0.01; % y方向采样间隔
u=min(xx)-0.5:du:max(xx)+0.5;
v=min(yy)-0.5:dv:max(yy)+0.5;
[X,Y]=meshgrid(u,v);
for i=1:mZ=real(local_density(i))*exp(-((X-xx(i))/sigma).^2-((Y-yy(i))/sigma).^2);Zs=Z+Zs;
end
mesh(X,Y,Zs)