3 files changed, 120 insertions, 0 deletions
diff --git a/figs/bands.fig/Makefile b/figs/bands.fig/Makefile
new file mode 100644
index 0000000..5a6a103
--- /dev/null
+++ b/figs/bands.fig/Makefile
@@ -0,0 +1,23 @@
+PROJECTNAME=bands
+
+all: $(PROJECTNAME).pdf
+
+$(PROJECTNAME).pdf:
+	python $(PROJECTNAME).py > $(PROJECTNAME).dat
+	gnuplot $(PROJECTNAME).gnuplot
+	pdflatex -jobname $(PROJECTNAME) $(PROJECTNAME)_plot.tex
+
+install: $(PROJECTNAME).pdf
+	cp $(PROJECTNAME).pdf $(INSTALLDIR)/
+
+
+clean-aux:
+	rm -f $(PROJECTNAME).aux
+	rm -f $(PROJECTNAME).log
+	rm -f $(PROJECTNAME).dat
+	rm -f $(PROJECTNAME)_plot.tex
+
+clean-pdf:
+	rm -f $(PROJECTNAME).pdf
+
+clean: clean-aux clean-pdf
diff --git a/figs/bands.fig/bands.gnuplot b/figs/bands.fig/bands.gnuplot
new file mode 100644
index 0000000..f1f84ff
--- /dev/null
+++ b/figs/bands.fig/bands.gnuplot
@@ -0,0 +1,42 @@
+#!/usr/bin/env gnuplot
+
+datafile="bands.dat"
+# output
+outfile="bands_plot.tex"
+set output outfile
+
+# terminal
+set term lua tikz size 10,6 standalone tightboundingbox
+
+# format axes
+set key off
+unset colorbox
+unset border
+unset xtics
+unset ytics
+unset ztics
+
+# colors
+## 4169E1 (pastel blue)
+## FF143C (bright red)
+## 32CD32 (bright green)
+## DAA520 (ochre)
+
+# angle camera
+set view 80,30
+
+# style for mesh lines
+set style line 1 linetype rgbcolor "#000000"
+
+# splot mode
+set pm3d hidden3d depthorder
+
+# set colors
+set palette defined (0 "#DAA520",0.33 "#FF143C", 0.66 "#32CD32",1 "#4169E1")
+
+# plot
+splot datafile using 1:2:3:4 with pm3d linestyle 1 , \
+ datafile using 5:6:7:8 with pm3d linestyle 1 , \
+ datafile using 9:10:11:12 with pm3d linestyle 1 , \
+ datafile using 13:14:15:16 with pm3d linestyle 1
+
diff --git a/figs/bands.fig/bands.py b/figs/bands.fig/bands.py
new file mode 100644
index 0000000..57d03b8
--- /dev/null
+++ b/figs/bands.fig/bands.py
@@ -0,0 +1,55 @@
+#!/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()