import pylab
import sys
import matplotlib.pyplot as plt

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

# The maximum temperature for the plot in K
Tmax = 400.0 

# The minimum temperature for the plot in K
Tmin = 250.0

# Calculated ideality factor
IF = 1.82080344438

# Voltages for which the calculations have been performed in V
voltages = [0.05,0.1,0.15,0.2,0.25,0.3]

# Name of output file containing the effective barrier associated with the thermionic emission process 
# from the Si(100) conduction band to Ag(100) in eV for each bias point
fname = 'schottkyeff'

#-------------------
# Physical constants
#-------------------

q = 1.60217662*10**-19
kB = 1.38064852*10**-23

#---------------------
# Plot initialization
#---------------------

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

#-----------------------------------------------------------
# Calculate the Schottky barrier and make the Arrhenius plot
#-----------------------------------------------------------

filenames = ['IvsT'+str(voltage)+'.dat' for voltage in voltages]

# Loop over the voltages
schottky = []
schottkyeff = []
for n in range(len(voltages)):

	# Get current and temperature
	IvsT = numpy.loadtxt(filenames[n], usecols=(0,1))
	x = 1./IvsT[1:,0] # x-values are 1/T
	y = (IvsT[1:,1])/(IvsT[1:,0]**2) # y-values are I/(T**2)
	
	# Fit a first degree polynomial to the values in the chosen temperature range
	m, b = numpy.polyfit(x[(x > 1./Tmax) & (x < 1./Tmin)], numpy.log(y[(x > 1./Tmax) & (x < 1./Tmin)]), 1)
	yfit = numpy.polyval([m,b],x)
	yfit2 = 10**(yfit/2.303) # Conversion from the natural log to log-10

	# For all the bias voltages considered calculate: 
	# 1) Schottky barrier
	# 2) Effective barrier associated with the thermionic emission process 
	#    from the Si(100) conduction band to Ag(100)
	schottky.append((voltages[n]/IF)-(kB/q)*m)
	schottkyeff.append(-(kB/q)*m)
	if schottkyeff[n] < 0:
		print()
		print("Bias = "+str(voltages[n])+" V")
		print("Schottky barrier =",schottky[n]*1000,"meV")
		print("No effective barrier")
 	elif schottkyeff[n] >= 0:
 		print()
		print("Bias = "+str(voltages[n])+" V")
		print("Schottky barrier =",schottky[n]*1000,"meV")
		print("Effective barrier =",schottkyeff[n]*1000,"meV")

	# Plot figure
	m = float(len(voltages))
	ax.semilogy(x*1000, y, label='V$_\mathrm{bias}$ = '+str(voltages[n])+' V',color=plt.cm.jet(n/m),linewidth=0,marker='o',markevery=10,fillstyle='none',markersize=4)
	ax.semilogy(x*1000, yfit2, color=plt.cm.jet(n/m),linestyle='-',linewidth=0.5)

	# Get the minumum and maximum y-values
	if n == 0:
		ymin = min(y[(x > 1./Tmax) & (x < 1./Tmin)])
	if n == len(voltages)-1:
		ymax = max(y[(x > 1./Tmax) & (x < 1./Tmin)])

# Write to file
out = open(fname+".dat","w")
for n in range(len(voltages)):
	out.write("%f %f\n"%(voltages[n], schottkyeff[n]))

ax.spines['top'].set_linewidth(1)
ax.spines['bottom'].set_linewidth(1)
ax.spines['right'].set_linewidth(1)
ax.spines['left'].set_linewidth(1)
ax.xaxis.set_tick_params(width=0.5)
ax.yaxis.set_tick_params(width=0.5)
ax.tick_params(direction='in', pad=10)

ax.tick_params(axis='both', which='major', labelsize=10)
ax.set_ylabel(r'$I T^{-2}$ (A K$^{-2}$)', fontsize=15)
ax.set_xlabel(r'1000 T$^\mathrm{-1}$ (K$^{-1}$)', fontsize=15)

ax.set_xlim(1000/Tmax,1000/Tmin)
ax.set_ylim(ymin-0.99*ymin,ymax+0.5*ymax)
ax.legend(loc="lower center",ncol=2,fontsize=10,numpoints=1,frameon=False)

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