import pylab as pl
import numpy as np

#--------------------------------------------------------------------#
# Script to calculate bulk resistivity using the MD-Landauer method  |
#--------------------------------------------------------------------#
#                        INPUT SECTION                               |
#--------------------------------------------------------------------#
cs_area = 8.65127469112**2
lengths_list = [0,12.235,16.313,20.391]
resistances_list = [1482.9,1535.1,1555.2,1571.90]
output_filename = 'resistivity.png'

#--------------------------------------------------------------------#
#                    END OF INPUT SECTION                            |
#--------------------------------------------------------------------#
#Fit the data and extract the 1D resistivity (Ohm/Angstrom)
(resistivity_1d,contact_resistance)=np.polyfit(lengths_list,resistances_list,1)
yp=np.polyval([resistivity_1d,contact_resistance],lengths_list)

#Convert the 1D resistivity from (Ohm/Angstrom) to (Ohm/Meter)
resistivity_1d = resistivity_1d*(10.0**10)
#Calculate the Bulk resistivity (Ohm*Meter)
resistivity_bulk = resistivity_1d * cs_area/((10.0**10)**2)

# Print report
print('+------------------------------------------------------------------------------+')
print('| Resistivity Report (T=300 K)                                                 |')
print('+------------------------------------------------------------------------------+')
print('| Contact resistance: %.2f Ohm' %(contact_resistance))
print('| 1D resistivity: %.2e Ohm/Meter' %(resistivity_1d))
print('| Bulk resistivity: %.2e Ohm*Meter' %(resistivity_bulk))
print('+------------------------------------------------------------------------------+')

# Plot the results
pl.figure()
ax = pl.subplot(111)
ax.tick_params(labelsize=14)
ax.set_xlabel(r"Length of the MD region ($\mathrm{\AA}$)",fontsize=20)
ax.set_ylabel(r"Resistance ($\Omega$)",fontsize=20)
ax.plot(lengths_list,resistances_list,marker='d',color='r',linestyle='None', label='300 K')
ax.plot(lengths_list,yp,linestyle='-',linewidth=1.5,color='r')
ax.text(0.35,0.25,r"$R_{\rm c}$ = %6.2f $\Omega$"%(contact_resistance),transform=ax.transAxes,fontsize=20)
ax.text(0.35,0.25-0.075,r"$\rho_\mathrm{1D}(T=300)$ = %.2e $\Omega / m$"%(resistivity_1d),transform=ax.transAxes,fontsize=20)
ax.text(0.35,0.25-2*0.075,r"$\rho_\mathrm{bulk}(T=300)$ = %.2e $\Omega \cdot m$"%(resistivity_bulk), transform=ax.transAxes,fontsize=20)
pl.tight_layout()
pl.savefig(output_filename,dpi=200)
pl.show()