import pylab
import sys
import matplotlib.pyplot as plt

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

# Filename of .nc file containing the device configuration and IVCurve analysis
fname = 'device_120.nc'

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

iv = nlread(fname, IVCurve)[0]
y = iv.currents().inUnitsOf(Ampere)
x = iv.biases().inUnitsOf(Volt)

# sort the list after bias
sort_ind = numpy.argsort(x)
x = numpy.sort(x)
y = y[sort_ind]

# Get the area perpendicular to the transport direction
config = nlread(fname,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
y2 = y/Area # Calculate current density

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

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

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.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax.tick_params(direction='in', pad=10)

ax.set_ylabel(r'$I$ (A $\mu$m$^{-2}$)', fontsize=15)
ax.set_xlabel(r'$V_{bias}$ (V)', fontsize=15)
ax.set_ylim(-max(y2),max(y2))
ax.set_xlim(float(min(x))-0.02,float(max(x))+0.02)

ax.plot([0,0],[-max(y2),max(y2)],color='black',linewidth=0.5,linestyle='--')
ax.plot([float(min(x))-0.02,float(max(x))+0.02],[0,0],color='black',linewidth=0.5,linestyle='--')

ax.plot(x, y2, marker='s',color='green',linewidth=1, markersize=4)
plt.tight_layout()
plt.savefig('IV.png',dpi=300)