clear all
close all

%Coordinates along kz axis
kx1=0;
ky1=0;
kz1=-0.5:0.001:0.5;
kz2=0:0.001:0.5;

%Coordinates along kx,ky at K and H points
k=-0.015:0.0001:0.015;
%k=-0.002:0.0001:0.002;
[kx,ky]=meshgrid(k);
kzK=0;
kzH=0.5;

%Band energies along kz
E1=Eband(kx1,ky1,kz1,1)
E2=Eband(kx1,ky1,kz1,2);
E3=Eband(kx1,ky1,kz1,3);
E4=Eband(kx1,ky1,kz1,4);
E5=Eband(kx1,ky1,kz2,1);
E6=Eband(kx1,ky1,kz2,2);
E7=Eband(kx1,ky1,kz2,3);
E8=Eband(kx1,ky1,kz2,4);


%Band energies at kz=0 (K point)
E2K=Eband(kx,ky,kzK,2);
E3K=Eband(kx,ky,kzK,3);

%Band energies at kz=0.5 (H point)
E2H=Eband(kx,ky,kzH,2);
E3H=Eband(kx,ky,kzH,3);

%Plotting of Band energies along kz
figure(1)
plot(kz1,E1(:),'LineWidth',2)
hold on
plot(kz1,E2(:),'LineWidth',2)
plot(kz1,E3(:),'LineWidth',2)
plot(kz1,E4(:),'LineWidth',2)
xlabel('kz (2\pi/c_0)')
ylabel('Energie (eV)')
fig = gcf;
set( findall(fig, '-property', 'fontsize'), 'fontsize', 18)

%Plotting detail of degenerate Band energies along kz
figure(2)
plot(kz2,1000*E5(:),'LineWidth',2)
hold on
plot(kz2,1000*E6(:),'LineWidth',2)
plot(kz2,1000*E7(:),'LineWidth',2)
plot(kz2,1000*E8(:),'LineWidth',2)
xlabel('kz (2\pi/c_0)')
ylim([-50 10])
ylabel('Energie (meV)')
fig = gcf;
set( findall(fig, '-property', 'fontsize'), 'fontsize', 18)

%Plotting of band energies at K point
figure(3)
surf(kx,ky,1000*E2K,'edgecolor', 'none');
hold on
[M1,c1]=contour3(kx,ky,1000*E2K,[-53.15,-53.15],'-b');
%[M2,c2]=contour3(kx,ky,1000*E2K,[-50,-50],'-r');
surf(kx,ky,1000*E3K,'edgecolor','none');
[M2,c2]=contour3(kx,ky,1000*E3K,[-40.6, -40.6],'-b');
c1.LineWidth=2;
c2.LineWidth=2;
colormap hot
colorbar
xlabel('kx (2\pi/c0)')
ylabel('ky (2\pi/c0)')
zlabel('Energie (meV)')

%Plotting of band energies at H point
figure(4)
surf(kx,ky,1000*E2H,'edgecolor', 'none');
hold on
[M1,c3]=contour3(kx,ky,E2H,[0,0],'-b');
surf(kx,ky,1000*E3H,'edgecolor', 'none');
[M1,c4]=contour3(kx,ky,E3H,[0,0],'-b');
colormap hot
colorbar
xlabel('kx (2\pi/c0)')
ylabel('ky (2\pi/c0)')
zlabel('Energie (meV)')

function E_N=Eband(kx,ky,kz,N)
  % Gets the eigenvalues of each hamiltonian matrix for each component of kx,ky,kz
  % Input: wave vectors kx,ky,ky normalized by 2*pi/c0, c0= 3.35A
  % supposes that kx and ky are generated by [kx,ky]=meshgrid(...)
  % N=1,2,3,4: each band index
  % Output: Matrix of dimensions Lx*Ly*Lz with the energies of the band N.
  % being Lx the length of kx taken as a vector 
  
  L=size(kx); % since kx is generated by meshgrid, kx is a matrix of dimensions Lx*Ly
  Lx=L(1);
  Ly=L(2);
  Lz=length(kz);
  
  %Initializing matrix
  E_N=zeros(L(1),L(2),Lz);
  
  % Calculating energies
  for kk=1:Lz
    for jj=1:Ly
        for ii=1:Lx
            dummyE=Hamiltonian(kx(ii,jj),ky(ii,jj),kz(kk));
            E_N(ii,jj,kk)= dummyE(N);
        end
    end
  end        
end 

function E_N=Hamiltonian(kx,ky,kz)
  % Hamiltonian of graphite considering hopping parameter gamma_3!=0
  % Input: 
  % kx,ky,kz: wave vectors normalized by 2*pi/c0, c0= 3.35A 
  % Output: E_N= Energies of the 4 bands. length(E_N)==4
  
  % Constants
  a0=2.46e-10;    % distance between carbon atoms in single graphite lattice
  e=1.6e-19;      % elementary charge
  hbar=1.055e-34; % Reduced Planck constant
  h=hbar*2*pi;    % Planck constant
  muB=5.789e-5;   % Bohr magneton
  gs = 2.5;       % value for g factor
 
  % Hopping parameters,
  g0= 3.37;
  g1= 0.363;
  g2= -0.0243;
  g3= 0.31; 
  g4= 0.07;
  g5= 0.05;
  Delta= -0.007; % It corresponds to gamma_6
  
  %Normalizing kx and ky
  c0=3.35e-10;
  kx=kx*2*pi/c0;
  ky=ky*2*pi/c0;
  
  GAMMA= 2*cos(kz*pi);
  sigma= sqrt(3).*a0.*sqrt(kx.^2 +ky.^2)./2;
  
  alpha=alphaf(kx,ky);
 
  e1= Delta +g1*GAMMA +0.5*g5*GAMMA.^2;
  e2= Delta -g1*GAMMA +0.5*g5*GAMMA.^2;
  e3= 0.5*g2*GAMMA.^2;
  
  H_13= 1/sqrt(2).* (-g0+g4.*GAMMA).*sigma.*exp(1j.*alpha);
  H_23= 1/sqrt(2).*(g0+g4.*GAMMA).*sigma.*exp(1j.*alpha);
  H_33= g3.*GAMMA.*sigma.*exp(1j.*alpha);
    
  H= [e1,0,H_13,conj(H_13); 0,e2,H_23,-conj(H_23); conj(H_13),conj(H_23),e3,H_33; H_13,-H_23,conj(H_33),e3];
 
  [E_N]= eig(H);
  E_N=sort(E_N);
  
end 

function a=alphaf(kx,ky)
   % Angle with respect to the ky axis
  if ky==0 & kx==0
      a=0;
  elseif kx>=0 & ky>=0 
      a=atan(-kx/ky); 
  elseif kx<=0 & ky>=0
      a=atan(-kx./ky);
  elseif kx<=0 & ky<=0
      a=pi+atan(-kx/ky);
  elseif kx>=0 & ky<=0
      a=-pi+atan(-kx/ky);
  end
end