% The Zhang lattice is a simple lattce that naturally has type-II Dirac points
% It was created by Yiqi Zhang's group in 2020
% This is the code for solving for the band structure of the Zhang lattice
% by using the plane-wave expansion method according to Eq. (3)
% zhangyiqi@xjtu.edu.cn

clc
clear all
close all

a=3;
pin=-5;
d1=0.25;
d2=0.75;

ym=(4+sqrt(3))*a;
xm=a;

xmax=15;
xN=2^10;
dx=2*xmax/xN;
x0=-xmax:dx:xmax-dx;

ymax=15;
yN=2^10;
dy=2*ymax/yN;
y0=-ymax:dy:ymax-dy;
[x,y]=meshgrid(x0,y0);

theta=-acos(1/(2*sqrt(5+2*sqrt(3))));
x1=cos(theta)*x-sin(theta)*y;
y1=sin(theta)*x+cos(theta)*y;
aL = sqrt((129+72*sqrt(3))/(20+8*sqrt(3)));
bL = (4+sqrt(3))/(2*sqrt(5+2*sqrt(3)));
cL = sqrt(5+2*sqrt(3));
x2=x1-aL*y1/bL;
y2=cL*y1/bL;

R=zeros(yN,xN);
for flagx=-5:5
    for flagy=-4:4        
            R0=pin*exp(-((x-flagx*xm).^2/d1^2+(y-flagy*ym).^2/d2^2))+...
                pin*exp(-((x-flagx*xm).^2/d1^2+(y-a-flagy*ym).^2/d2^2))+...
                pin*exp(-((x-flagx*xm).^2/d1^2+(y+a-flagy*ym).^2/d2^2))+...
                pin*exp(-((x-a/2-flagx*xm).^2/d1^2+(y-a-sqrt(3)*a/2-flagy*ym).^2/d2^2))+...
                pin*exp(-((x-a/2-flagx*xm).^2/d1^2+(y-2*a-sqrt(3)*a/2-flagy*ym).^2/d2^2))+...
                pin*exp(-((x-a/2-flagx*xm).^2/d1^2+(y-3*a-sqrt(3)*a/2-flagy*ym).^2/d2^2));
       
        R=R+R0;
    end
end

Dx = cL*a;
Dy = cL*a;
xmax=Dx/2;
xN=2^10;
dx=2*xmax/xN;
x0=-xmax:dx:xmax-dx;
ymax=Dy/2;
yN=2^10;
dy=2*ymax/yN;
y0=-ymax:dy:ymax-dy;
[xin,yin]=meshgrid(x0,y0);
Rin=griddata(x2,y2,R,xin,yin);

Kx=2*pi/Dx;
Ky=2*pi/Dy;
pN=15;
qN=15;
flagx=1;
for p=-pN:pN
    flagy=1;
    for q=-qN:qN
        vn(flagy,flagx)=sum(sum(Rin.*exp(-1i*Kx*p*xin-1i*Ky*q*yin)))*dx*dy/(Dx*Dy);
        flagy=flagy+1;
    end
    flagx=flagx+1;
end

kN=200;
dkx=Kx/kN;
dky=Ky/kN;
Sx=pN;
Sy=qN;
Nlx=2*Sx+1;
Nly=2*Sy+1;
MV=zeros(Nlx*Nly);
for flag1=-pN:pN
    Mi=zeros(Nly);
    for flag2=-qN:qN
        Mi=Mi+...
            diag(vn(flag2+qN+1,flag1+pN+1)*ones(1,Nly-abs(flag2)),-flag2);
    end
    if flag1<=0
        for ind=0:2*Sx+flag1
            MV(ind*Nly+[1:Nly],(-flag1+ind)*Nly+[1:Nly])=Mi;
        end
    else
        for ind=0:2*Sx-flag1
            MV((flag1+ind)*Nly+[1:Nly],ind*Nly+[1:Nly])=Mi;
        end
    end
end
kx=-4*Kx:dkx:4*Kx;
ky=-4*Ky:dky:4*Ky;
eigval1=zeros(length(kx));
eigval2=zeros(length(kx));
eigval3=zeros(length(kx));
for indx=1:length(kx)
    k1=kx(indx);
    for indy=1:length(kx)
        fprintf('indx=%2d, indy=%2d \n',indx,indy);
        k2=ky(indy);
        
        Mtemp=zeros(Nly,Nlx);
        for m=-Sx:Sx
            for n=-Sy:Sy
                Mtemp(n+Sy+1,m+Sx+1)=(1/2+1/2*(aL/bL)*(aL/bL))*(Kx*m+k1)^2+...
                    (1/(2*(bL/cL)*(bL/cL)))*(Ky*n+k2)^2-(aL*cL/(bL*bL))*(Kx*m+k1)*(Ky*n+k2);
            end
        end
        G=reshape(Mtemp,Nlx*Nly,1);
        M=MV+diag(reshape(Mtemp,Nlx*Nly,1));
        eigenvalues=sort(real(eig(M)));
        
        eigval1(indx,indy)=eigenvalues(1);
        eigval2(indx,indy)=eigenvalues(2);
        eigval3(indx,indy)=eigenvalues(3);
    end
end

[kx,ky]=meshgrid(kx,ky);
KX=kx;
KY=-aL/bL*kx+cL/bL*ky;
kx1=cos(-theta)*KX-sin(-theta)*KY;
ky1=sin(-theta)*KX+cos(-theta)*KY;

figure('Units','pixels','Position',[128 128 1000 1000],...
    'PaperUnits','inches','PaperPositionMode','manual',...
    'PaperPosition',[0 0 10 10]);
hold on
surf(kx1,ky1,-eigval1)
surf(kx1,ky1,-eigval2)
surf(kx1,ky1,-eigval3)
shading interp
camlight head
hold off
axis equal
axis([-pi/a-0.01 pi/a+0.01 -pi/2/a-0.005 pi/2/a+0.005 0.2 0.62])
daspect([1 1 0.45/8.5*a])
view(10,25)
box on
