import pylab
import sys
import matplotlib.pyplot as plt

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

# Name of the input file containing the device configuration and I-V data
filename = 'device_120.hdf5'

# Temperature at which calculations have been run (in Kelvin)
T = 300

# Number of forward bias points to be included in the fit
numpoints = 2

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

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

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

iv = nlread(filename, IVCurve)[0]
y = numpy.array(iv.currents())
x = numpy.array(iv.biases())

# Get the area perpendicular to the transport direction
config = nlread(filename,DeviceConfiguration)[1]
a = config.bravaisLattice().primitiveVectors()[0][0].inUnitsOf(Ang)*10**(-4)
b = config.bravaisLattice().primitiveVectors()[1][1].inUnitsOf(Ang)*10**(-4)
Area = a*b

#---------------------------
# Calculate ideality factor  
# --------------------------

# Exclude the zero bias point
y2 = y[abs(y) > 1e-20]
x2 = x[x != 0]

# Get y-values for fit
logy = numpy.log(abs(y2)/(1-numpy.exp((-q*abs(x2))/(kB*T))))

# Consider only positive bias
xpos = x2[x2 > 0] 
logypos = logy[x2 > 0]

# Fit to a first degree polynomial and print out the ideality factor
(m_n,b_n) = numpy.polyfit(xpos[0:numpoints],logypos[0:numpoints],1) 
yp = numpy.polyval([m_n,b_n],x)
idfact = q/(kB*T*m_n)

print("+------------------------------------------------------------------------------+")
print()
print("Forward bias points considered in the fit:")
for i in numpy.arange(0,numpoints):
    print(xpos[i], "V")
print()
print("Ideality factor n =",idfact)
print()
print("-------------------------------------------------------------------------------+")

#-------------
# Plot figure 
#-------------

plt.figure(figsize=(5,5))
ax = pylab.subplot(111)

yp2 = 10**(yp/2.303) # Convert from the natural log to log-10

plt.semilogy(x2, abs(y2)/(1-numpy.exp((-q*abs(x2))/(kB*T)))/Area, marker='s',color='green',linewidth=0,markersize=6)
plt.semilogy(x[x >= 0], yp2[x >= 0]/Area,color='green',linewidth=0.5,linestyle='-')
plt.semilogy([0,0],[10**-16,5*10**-7],color='black',linewidth=0.5)

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

ax.set_ylabel(r'$I/(1-e^{-qV_{bias}/k_BT}$) (A $\mu$m$^{-2}$)', fontsize=15)
ax.set_xlabel(r'$V_{bias}$ (V)', fontsize=15)
ax.set_ylim(min(yp2[x >= 0]/Area)-min(yp2[x >= 0])/(2*Area),max(yp2[x >= 0]/Area))
ax.set_xlim(0,max(x[x >= 0])+0.02)
plt.tight_layout()
plt.savefig('IV-n-log.png',dpi=300)