#!/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, ax = plt.subplots(nrows=2, ncols=1, figsize=(3.69,4.5))
#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_FS = np.array(
[[119.5, 0.262, 0.020, 0.999, 0.999]]
)

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

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[:,0], LO_475[:,3], yerr=LO_475[:,4], marker='o', color='b', linestyle='--', fillstyle='none' )
plt.errorbar( LO_1100[:,0], LO_1100[:,1], yerr=LO_1100[:,2], marker='o', color='r', label='AIMD 1100 K' )
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[:,0], LO_1100[:,3], yerr=LO_1100[:,4], marker='o', color='r', linestyle='--', fillstyle='none' )
plt.plot( Guo[:,0], Guo[:,1], marker='o', color='g', label='Hu et al.$\cite{hu_vibrational_2018}$' )

plt.annotate('(a)', xy=(110, 0.04))

plt.legend(loc='upper left', numpoints=1, frameon=False, fontsize=10)
plt.xlim(0,130.)
plt.ylim(0,0.3)
plt.tick_params(labelbottom='off', length=6, width=1)
plt.ylabel('Sticking probability')

plt.subplot(2,1,2)

plt.scatter( NTPD_475_E, np.array(NTPD_475_S)*1.5, marker='s', color='b', facecolor='none', linewidths=0.5 )
plt.scatter( Hdes_475_E, np.array(Hdes_475_S)*1.5, marker='d', color='b', facecolor='none', linewidths=0.5 )
plt.scatter( Ndes_1100_E, np.array(Ndes_1100_S)*1.5, marker='s', color='r', facecolor='none', linewidths=0.5 )
plt.scatter( Hdes_1100_E, np.array(Hdes_1100_S)*1.5, marker='d', color='r', facecolor='none', linewidths=0.5 )
plt.plot(x_axis, np.array(fit_475_S)*1.5, color='b', label=r'Exp. 475 K $\times$ 1.5', linestyle='--')
plt.plot(x_axis, np.array(fit_1100_S)*1.5, color='r', label=r'Exp. 1100 K $\times$ 1.5', linestyle='--')
plt.errorbar( LO_475[:,0], LO_475[:,1], yerr=LO_475[:,2], marker='o', label='AIMD 475 K', color='b' )
plt.errorbar( LO_1100[:,0], LO_1100[:,1], yerr=LO_1100[:,2], marker='o', label='AIMD 1100 K', color='r' )

plt.plot( [70., 74.], [0.046, 0.046], c='k' )
plt.annotate('4.0', xy=(79., 0.046), fontsize=8)
plt.plot( [87.8, 92.], [0.092, 0.092], c='k' )
plt.annotate('4.2', xy=(97., 0.092), fontsize=8)
plt.plot( [102.9, 110.], [0.146, 0.146], c='k' )
plt.annotate('7.1', xy=(115., 0.140), fontsize=8)
plt.plot( [118.5, 119.5], [0.172, 0.172], c='k' )
plt.annotate('1.0', xy=(121.5, 0.166), fontsize=8)

plt.plot( [102.9, 106.4], [0.132, 0.132], c='k' )
plt.annotate('3.5', xy=(90., 0.132), fontsize=8)
plt.plot( [119.5, 126.5], [0.194, 0.194], c='k' )
plt.annotate('7.0', xy=(107., 0.194), fontsize=8)

plt.annotate('(b)', xy=(110, 0.04))

plt.legend(loc='upper left', numpoints=1, frameon=False, fontsize=10)
plt.xlim(0,130.)
plt.ylim(0,0.249)
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.pgf')
#plt.show()