#!/usr/bin/env python
import matplotlib as mpl
mpl.use('pgf')
import matplotlib.pyplot as plt
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, Reaction + Trapping, error
LO_1100 = np.array(
[[74.0, 0.046, 0.009, 0.196, 0.018],
[87.8, 0.092, 0.013, 0.192, 0.018],
[102.9, 0.146, 0.016, 0.220, 0.019],
[119.5, 0.172, 0.017, 0.224, 0.019]]
)

LO_1100_v0 = np.array(
[[74.0, 0.04583333333333333, 0.009545138257346085, 0.19583333333333333, 0.018113222546485716],
[87.8, 0.08425720620842572, 0.01307982734569903, 0.18181818181818182, 0.018161649769418916],
[102.9, 0.13043478260869565, 0.01611042771695256, 0.20594965675057209, 0.01934478185394897],
[119.5, 0.16120906801007556, 0.018455522434108036, 0.22418136020151133, 0.02093074089003594]]
)

LO_1100_FS = np.array(
[[119.5, 0.262, 0.020, 0.999, 0.999]]
)

LO_1100_FS_v0 = np.array(
[[119.5, 0.251256281407, 0.0217412006379, 0.27135678392, 0.0222887785127]]
)

LO_475 = np.array(
[[102.9, 0.132, 0.011, 0.209, 0.013],
[119.5, 0.194, 0.013, 0.244, 0.014]]
)

LO_475_v0 = np.array(
[[102.9, 0.12060889929742388, 0.011144273831711437, 0.1955503512880562, 0.013572185852126157],
[119.5, 0.18523153942428036, 0.013743621400307245, 0.23779724655819776, 0.015061397751776505]]
)

Guo = np.array(
[[38.6, 0.01, 0.05],
[57.9, 0.03, 0.05],
[77.2, 0.08, 0.09],
[96.5, 0.14, 0.14],
[115.8,0.24, 0.24],
[135.1,0.33, 0.33],
[154.4,0.41, 0.41],
[173.7,0.48, 0.48]]
)

NTPD_475_E = [4.9,31.2,54.0,44.7,56.7,80.1,94.6,106.1,117.2]				# square
NTPD_475_S = [0.018,0.023,0.034,0.023,0.025,0.056,0.063,0.086,0.106]
Hdes_475_E = [4.9,31.4,44.8,56.8]							# diamond
Hdes_475_S = [0.019,0.026,0.022,0.024]
Hdes_1100_E = [4.9,31.4,39.6,56.7]							# triangle
Hdes_1100_S = [0.002,0.005,0.005,0.013]
Ndes_1100_E = [31.3,39.6,44.8,54.1,56.7,67.6,80.2,92.8,94.6,106.2,117.3,125.4,121.6]	# circle
Ndes_1100_S = [0.004,0.006,0.008,0.016,0.013,0.024,0.052,0.064,0.070,0.088,0.111,0.118,0.129]
E_475 = NTPD_475_E + Hdes_475_E
S_475 = NTPD_475_S + Hdes_475_S
E_1100 = Hdes_1100_E + Ndes_1100_E
S_1100 = Hdes_1100_S + Ndes_1100_S

#plt.subplot(2,1,1)

x_axis = np.linspace(0.,150.,150)
#plt.scatter( NTPD_475_E, NTPD_475_S, marker='s', color='b', facecolor='none', linewidths=0.5 )
#plt.scatter( Hdes_475_E, Hdes_475_S, marker='d', color='b', facecolor='none', linewidths=0.5 )

#plt.scatter( Ndes_1100_E, Ndes_1100_S, marker='s', color='r', facecolor='none', linewidths=0.5 )
#plt.scatter( Hdes_1100_E, Hdes_1100_S, marker='d', color='r', facecolor='none', linewidths=0.5 )

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', linestyle='--')

fit_1100_S = []
A = 0.4
E0 = 150.
W = 80.
for i in range(len(x_axis)):
	fit_1100_S.append( A/2. * ( 1 + erf((x_axis[i] - E0)/W) ) )
#plt.plot(x_axis, fit_1100_S, color='r', label='Exp. 1100 K', linestyle='--')

plt.errorbar( LO_475[:,0], LO_475[:,1], yerr=LO_475[:,2], marker='o', color='b', label='AIMD 475 K' )
plt.errorbar( LO_475_v0[:,0], LO_475_v0[:,1], yerr=LO_475_v0[:,2], marker='o', markerfacecolor='none', markeredgecolor='b', linestyle='--', color='b', label='' )
plt.errorbar( LO_1100[:,0], LO_1100[:,1], yerr=LO_1100[:,2], marker='o', color='r', label='AIMD 1100 K' )
plt.errorbar( LO_1100_v0[:,0], LO_1100_v0[:,1], yerr=LO_1100_v0[:,2], marker='o', markerfacecolor='none', markeredgecolor='r', linestyle='--', color='r', label='' )
plt.errorbar( LO_1100_FS[:,0], LO_1100_FS[:,1], yerr=LO_1100_FS[:,2], fmt='o', color='k', label='AIMD 1100 K (FS)' )
plt.errorbar( LO_1100_FS_v0[:,0], LO_1100_FS_v0[:,1], yerr=LO_1100_FS_v0[:,2], marker='o', markerfacecolor='none', markeredgecolor='k', linestyle='--', color='k', label='' )

plt.plot( Guo[:,0], Guo[:,1], marker='o', markerfacecolor='none', markeredgecolor='g', linestyle='--', color='g', label=r'Hu et al.$\textrm{\cite{hu_vibrational_2018}}$' )
#plt.plot( Guo[:,0], Guo[:,1], marker='o', markerfacecolor='none', markeredgecolor='g', linestyle='--', color='g', label=r'Hu et al.' )

plt.legend(loc='upper left', numpoints=1, frameon=False, fontsize=10)
plt.xlim(0,130.)
plt.ylim(0,0.3)
plt.tick_params(length=6, width=1)
plt.ylabel('Sticking probability')
plt.xlabel('Incidence energy (kJ/mol)')

plt.tight_layout()
plt.subplots_adjust(wspace=0, hspace=0)
plt.savefig('reactionprobability_v0.pgf')
#plt.show()