# 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("PBE_MSD_Coefficients.dat")
ds2<-read.table("SCAN_MSD_Coefficients.dat")

#Bin Size
bins1=6
bins2=4
#Range
xmin=0
xmax=4
ymin=0
ymax=13

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

hist(ds1[[1]],
main="",
xlim=c(xmin,xmax),
ylim=c(ymin,ymax),
xlab=expression(paste("D (10" ^-5*"cm"^2*"/s)")),
ylab="N",
col="#1DC9FF",
breaks=bins1,
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(paste("D (10" ^-5*"cm"^2*"/s)")),
ylab="N",
col="#FF0000",
breaks=bins2,
cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5, add=T)

pdf("MSD_Histogram.pdf", width=7, height=5)

hist(ds1[[1]],
main="",
xlim=c(xmin,xmax),
ylim=c(ymin,ymax),
xlab=expression(paste("D (10" ^-5*"cm"^2*"/s)")),
ylab="N",
col="#1DC9FF",
breaks=bins1,
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(paste("D (10" ^-5*"cm"^2*"/s)")),
ylab="N",
col="#FF0000",
breaks=bins2,
cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5, add=T)