# MakeFirgure.R: generate Histogram plot from epsilon infinity
# use: on the linux command line:
# $ Rscript MakeFirgure.R

# eps_inf*.dat: electric permitivity calculated from a_iso
# read data into a data frame
ds1<-read.table("PBE400.dat")
ds2<-read.table("SCAN330.dat")

# Number of breaks
brk1=7
brk2=7
#Range
xmin=1.76
xmax=1.90
ymin=5
ymax=280

setEPS()
postscript("eps_inf.eps", width=7, height=5)

hist(ds1[[1]],
main="",
xlim=c(xmin,xmax),
ylim=c(ymin,ymax),
xlab=expression(epsilon[infinity]),
ylab="N",
col="#1DC9FF",
breaks=brk1,
cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)
hist(ds2[[1]],
main="",
xlim=c(xmin,xmax),
ylim=c(ymin,ymax),
xlab=expression(epsilon[infinity]),
ylab="N",
col="#FF0000",
breaks=brk2,
cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5, add=T)
segments(1.77,ymin,1.77,ymax-35,col="black",lwd=4)
segments(1.813,ymin+125,1.813,ymax-35,col="#FF0000",lwd=4)
segments(1.866,ymax-55,1.866,ymax-35,col="#1DC9FF",lwd=4)
text(x=1.774, y=ymax-05, labels = "Exp. (1.77)",cex=1.5)
text(x=1.82, y=ymax-05, labels = "SCAN (1.81)",cex=1.5)
text(x=1.87, y=ymax-05, labels = "PBE (1.87)",cex=1.5)