clear;
hbar2_d_2m=1;
GX=[-1 -1 -1  0 0 0  1 1 1]; %in Units of 2.*pi/a_lat
GY=[-1  0  1 -1 0 1 -1 0 1]; %in Units of 2.*pi/a_lat
V= [ 0.00  0.2  0.0  0.2 0 0.2  0.0 0.2 0.0];
VMatrix=zeros(length(GX));
for g=1:length(GX)
    for q=1:length(GX)
        g_prime=find((GX==(GX(g)+GX(q))).*(GY==(GY(g)+GY(q))));
        VMatrix(g,g_prime)=VMatrix(g,g_prime)+V(q);
    end
end
[kx,ky]=meshgrid(-0.5:0.02:0.5); %in Units of 2.*pi/a_lat;
for a=1:size(kx,1);
    for b=1:size(kx,2)
        
    H=hbar2_d_2m.*diag((GX+kx(a,b)).^2+(GY+ky(a,b)).^2)+VMatrix;
    [EV,EE]=eig(H);
    EVal(a,b,:)=diag(EE); 
    end

end
contour(kx,ky,EVal(:,:,1),30)