#!/usr/bin/env python

## compute the bands for bilayer graphene (no g4 or delta)

from math import *
import cmath
import numpy

g1=0.1
g3=0.33*g1

def Omega(x,y):
    return(1+2*cmath.exp(-3j/2*x)*cos(sqrt(3)/2*y))

# Hamiltonian
def H(x,y):
    return(numpy.array(\
                       [[0,g1,0,Omega(x,y).conjugate()],\
                        [g1,0,Omega(x,y),0],\
                        [0,Omega(x,y).conjugate(),0,g3*Omega(x,y)*cmath.exp(3j*x)],\
                        [Omega(x,y),0,g3*(Omega(x,y).conjugate())*cmath.exp(-3j*x),0]]\
    ))

# eigenvalues
def eigsH(x,y):
    return(numpy.linalg.eigvals(H(x,y)))


# resolution
nrpoints=55
# xrange
xmin,xmax=0,4*pi/3
# yrange
ymin,ymax=-2*pi/sqrt(3),2*pi/sqrt(3)

# sample points
x=numpy.linspace(xmin,xmax,nrpoints)
y=numpy.linspace(ymin,ymax,nrpoints)
x,y=numpy.meshgrid(x,y)

# data points
z=numpy.zeros((4,nrpoints,nrpoints))
for i in range(0,nrpoints):
    for j in range(0,nrpoints):
        eigs=numpy.sort(numpy.real(eigsH(x[i,j],y[i,j])))
        for k in range(0,4):
            z[k,i,j]=(eigs[k])

# output
for i in range(0,nrpoints):
    for j in range(0,nrpoints):
        for k in range(0,4):
            print("%10f %10f %10f %d   " % (x[i,j],y[i,j],z[k,i,j],k),end='')
        print()
    print()