clear;
hbar2_d_2m=1;
G=[-5 -4 -3 -2 -1 0 1 2 3 4 5]; %in Units of 2.*pi/a_lat
V=[0 0.00 0 0.1 -0.1 0 -0.1 0.1 0 0.00 0];
VMatrix=zeros(length(G));
for g=1:length(G)
    for q=1:length(G)
        g_prime=find(G==(G(g)+G(q)));
        VMatrix(g,g_prime)=VMatrix(g,g_prime)+V(q);
    end
end
%imagesc(VMatrix)
kreal=[0:0.01:0.5]; %in Units of 2.*pi/a_lat;
kappa=[0:0.01:0.2].*i;
k=[fliplr(kappa),kreal,0.5+kappa];
for a=1:length(k);
    H=hbar2_d_2m.*diag((G+k(a)).^2)+VMatrix;
    [EV,EE]=eig(H);
    EVal(a,:)=sortrows(diag(EE)); 
    
end



figure(1)
subplot(1,2,1)
plot(real(EVal),'.')
hold on
plot([length(kappa) length(kappa)],[-0.5 3])
plot([length(kappa)+length(kreal) length(kappa)+length(kreal)],[-0.5 3])
ylim([-0.50 3])
hold off
subplot(1,2,2)
plot(abs(imag(EVal)),'.')
hold on
plot([length(kappa) length(kappa)],[-0.5 3])
plot([length(kappa)+length(kreal) length(kappa)+length(kreal)],[-0.5 3])

ylim([-0.50 3])
hold off