#!/usr/bin/env python
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
from math import erf

mpl.rc("text", usetex=True)
fig = plt.figure(figsize=(3.69,3.))
#plt.rcParams.update({'font.size': 10})

#Ei, Reaction, error
NN_ds = np.array(
[[0.94,0.002491,0.001759],#still running
[1.18,0.031950,0.001244],
[1.29,0.046900,0.001495],
[1.55,0.127350,0.002357],
[1.80,0.195260,0.002803],
[2.12,0.267403,0.003130],
[2.56,0.343756,0.003359]]
)

NN_fs = np.array(
[[0.94, 0.001000, 0.000224],
[1.18, 0.025300, 0.001110],
[1.29, 0.039050, 0.001370],
[1.55, 0.113900, 0.002246],
[1.80, 0.182650, 0.002732],
[2.12, 0.255388, 0.003084],
[2.56, 0.335300, 0.003338]]
)

NN = np.array(
[[0.94, 0.001150, 0.000240],
[1.18, 0.019303, 0.000973],
[1.29, 0.033152, 0.001266],
[1.55, 0.105661, 0.002174],
[1.80, 0.171730, 0.002668],
[2.12, 0.245175, 0.003046],
[2.56, 0.335325, 0.003345]]
)

NN_v2 = np.array(
[[0.28, 0.000888, 0.000096],
[0.32, 0.002789, 0.000169],
[0.45, 0.037910, 0.000616],
[0.52, 0.074104, 0.000847],
[0.78, 0.190009, 0.001277],
[0.97, 0.249733, 0.001444],
#[[0.41, 0.027999, 0.002123],
#[0.61, 0.136348, 0.003704],
#[0.82, 0.223456, 0.004148], 
#[0.94, 0.260376, 0.003103],
[1.18, 0.313047, 0.003279],
[1.29, 0.339317, 0.003348],
[1.55, 0.392707, 0.003454],
[1.80, 0.429594, 0.003502],
[2.12, 0.466750, 0.003533],
[2.56, 0.511784, 0.003547]]
)

Shirhatti = np.array(
[[0.94, 0.000006, 0.000005],
[1.18, 0.00003, 0.00001],
[1.29, 0.00012, 0.00007],
[1.55, 0.0012, 0.0001],
[1.80, 0.0045, 0.0004],
[2.12, 0.0082, 0.0008],
[2.56, 0.021, 0.007]]
)

SRP32vdW = np.array(
[[1.29, 0.162, 0.017],
[1.80, 0.266, 0.020],
[2.56, 0.382, 0.022]]
)

#print( 'Distorted surface/mobile surface:', NN_ds / NN )
#print( 'Frozen surface/mobile surface:', NN_fs / NN )

#fit_475_S = []
#A = 0.4
#E0 = 160.
#W = 80.
#for i in range(len(x_axis)):
#	fit_475_S.append( A/2. * ( 1 + erf((x_axis[i] - E0)/W) ) + 0.018 )
#plt.plot(x_axis, fit_475_S, color='b', label='Exp. 475 K')

#plt.errorbar( Shirhatti[:,0]*96.48530749926, Shirhatti[:,1]*6., yerr=Shirhatti[:,2], marker='s', label='Exp. (x6)', color='g', capsize=4. )
#plt.errorbar( SRP32vdW[:,0]*96.48530749926, SRP32vdW[:,1], yerr=SRP32vdW[:,2], marker='d', label='SRP32-vdW-DF1', color='k', capsize=4. )
#plt.errorbar( NN_fs[:,0]*96.48530749926, NN_fs[:,1], yerr=NN_fs[:,2], marker='o', label='MS.RPBE (FS)', color='b', capsize=4. )
#plt.errorbar( NN_ds[:,0]*96.48530749926, NN_ds[:,1], yerr=NN_ds[:,2], marker='o', label='MS.RPBE (DS)', color='dimgray', capsize=4. )
plt.errorbar( NN[:,0]*96.48530749926, NN[:,1], yerr=NN[:,2], marker='o', label='Molecular beam', color='r', capsize=4. )
plt.errorbar( NN_v2[:,0]*96.48530749926, NN_v2[:,1], yerr=NN_v2[:,2], marker='o', label=r'$\nu=2,J=1$', color='orange', capsize=4. )

dErovib = 66.9822596142

dx = 122.
print( dx / dErovib )
i = 6
plt.annotate('{:4.1f}'.format(dx), xy=(NN[i,0]*96.48530749926 - 0.5 * dx, NN[i,1]+0.01), size=10)
plt.plot( [NN[i,0]*96.48530749926 - dx, NN[i,0]*96.48530749926], [NN[i,1],NN[i,1]], zorder=10, c='k' )

dx = 113.
print( dx / dErovib )
i = 5
plt.annotate('{:4.1f}'.format(dx), xy=(NN[i,0]*96.48530749926 - 0.5 * dx, NN[i,1]+0.01), size=10)
plt.plot( [NN[i,0]*96.48530749926 - dx, NN[i,0]*96.48530749926], [NN[i,1],NN[i,1]], zorder=10, c='k' )

dx = 103.
print( dx / dErovib )
i = 4
plt.annotate('{:4.1f}'.format(dx), xy=(NN[i,0]*96.48530749926 - 0.5 * dx, NN[i,1]+0.01), size=10)
plt.plot( [NN[i,0]*96.48530749926 - dx, NN[i,0]*96.48530749926], [NN[i,1],NN[i,1]], zorder=10, c='k' )

dx = 93.
print( dx / dErovib )
i = 3
plt.annotate('{:4.1f}'.format(dx), xy=(NN[i,0]*96.48530749926 - 0.5 * dx, NN[i,1]+0.01), size=10)
plt.plot( [NN[i,0]*96.48530749926 - dx, NN[i,0]*96.48530749926], [NN[i,1],NN[i,1]], zorder=10, c='k' )

dx = 83.
print( dx / dErovib )
i = 2
plt.annotate('{:4.1f}'.format(dx), xy=(NN[i,0]*96.48530749926 - 0.5 * dx, NN[i,1]+0.01), size=10)
plt.plot( [NN[i,0]*96.48530749926 - dx, NN[i,0]*96.48530749926], [NN[i,1],NN[i,1]], zorder=10, c='k' )

plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(0., 270.)
plt.ylim(0.0, 0.6)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True)
plt.ylabel('Reaction probability')
plt.xlabel('Incidence energy (kJ/mol)')

plt.tight_layout()
plt.savefig('reactionprobability_vibrational.pdf')
#plt.show()