import pylab
import sys
import matplotlib.pyplot as plt

#------------------------
# User defined parameters
#------------------------

# The .dat output file from arrhenius.py
fname = 'schottkyeff'

# The .dat output file containing the barriers associated with the thermionic emission 
# process in eV for each bias point from spectral_current.py
fname2 = 'phi_F'

# The .dat output file containing the energy of maximum spectral current in eV for each 
# bias point from spectral_current.py
fname3 = 'Emax'

# The ideality factor
n = 1.82080344438

# The Schottky barrier calculated with pldos_hdp.py at zero bias
phipot = 0.606

#-----------
# Load files
#-----------

data = numpy.loadtxt(fname+'.dat', usecols=(0,1))
V = data[:,0]
Phi_F = data[:,1]

data = numpy.loadtxt(fname2+'.dat', usecols=(0,1))
V2 = data[:,0]
Phieff = data[:,1]
Vmax = max(V2)

data = numpy.loadtxt(fname3+'.dat', usecols=(0,1))
V3 = data[:,0]
Emax = data[:,1]

#------------
# Plot Figure
#------------

fig = plt.figure(figsize=(4,3))
ax = fig.add_subplot(111)

m = float(len(V))
ax.plot([0,(Vmax+0.02)],[phipot,phipot-(Vmax+0.02)/n],color='k',linewidth=.5,linestyle='--')

ax.plot(V2,Phi_F,color='k',fillstyle='none',linestyle='-',linewidth=.5)
for i in range(len(V)):
     ax.scatter(V[i],Phi_F[i],marker='s',facecolor=plt.cm.jet(i/m),zorder=10,s=20,edgecolor='k',linewidth=.5)

for i in range(len(V2)):
     ax.scatter(V2[i],Phieff[i],facecolor=plt.cm.jet(i/m),marker='^',zorder=10,s=20,edgecolor='k',linewidth=.5)

ax.plot(V3,Emax-V3/2,color='k',linewidth=.5)
for i in range(len(V3)):
     ax.scatter(V3[i],Emax[i]-V3[i]/2,facecolor=plt.cm.jet(i/m),zorder=10,s=20,edgecolor='k',linewidth=.5)

ax.set_xlim(0,max(V)+0.02)
ax.set_xlabel(r'$V_{bias}$ (V)', fontsize=10)
ax.set_ylabel('$E$-$\mu_R$ (eV)', fontsize=10)

plt.tight_layout()
plt.savefig('barrier_compare.png',dpi=300)