clear
% Anzahl der k-Punkte pro Richtung
nk=121;                           
% Gitterkonstante
lc=1;             
% Vektoren zum nächsten Nachbarn
delta_1=lc.*[sqrt(3)/2, 0.5];                                               
delta_2=lc.*[-sqrt(3)/2,0.5];            
delta_3=lc.*[0, -1];
% Hybridisierung; "Hopping" parameter
t=1;                                 
% Abstand des K-Punktes vom Gamma Punkt im k-Raum
K=(4/3/sqrt(3).*pi/lc);                                                      
% Erzeugen eines Netzes von äquidistanten Punkten im k-Raum
[kx,ky]=meshgrid(-1.5.*K:(3.*K/(nk-1)):1.5.*K);                             

% Für jeden k-Punkt wird jetzt eine 2x2 Hamiltonmatrix berechnet und diagonalisiert:
for a=1:size(kx,1)
    for b=1:size(kx,2)
        % Skalarprodukt von k und delta_1
        k_delta_1=kx(a,b)*delta_1(1)+ky(a,b)*delta_1(2);
        % Skalarprodukt von k und delta_2
        k_delta_2=kx(a,b)*delta_2(1)+ky(a,b)*delta_2(2);     
        % Skalarprodukt von k und delta_3
        k_delta_3=kx(a,b)*delta_3(1)+ky(a,b)*delta_3(2);                    
        Delta=exp(i.*k_delta_1)+exp(i.*k_delta_2)+exp(i.*k_delta_3);   
        % Hamilton-Matrix (2x2) pro k-Punkt
        H=([0                       t.*Delta;                               
            conj(t).*conj(Delta)    0           ]);
        % diagonalisieren V: Eigenvektoren, E: Eigenwerte
        [V,E]=eig(H);  
        % "Wegschreiben der Energien: E1 erstes Band, E2 zweites Band
        E1(a,b)=E(1,1);                                                     
        E2(a,b)=E(2,2);
    end
end
% zeichne erstes Band
mesh(kx,ky,E1)                                                              
hold on
% zeichne zweites Band
mesh(kx,ky,E2)  
% zeichne Brillouin Zone
phi=[0:6].*pi/3;
plot3(K.*cos(phi),K.*sin(phi),0.*phi,'k')                                   
axis equal
hold off                                                                    
% fertig