#!/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,4.5))
#plt.rcParams.update({'font.size': 10})

#2 indicates that trapping is included as well
LO_SRP32 = np.array(
[[119.5, 0.205, 0.0285]]
)

#Ei, Reaction, error, Trapping, error
LO_SRP100 = np.array(
[[74.0, 0.048, 0.010, 0.150, 0.016],
[87.8, 0.092, 0.013, 0.100, 0.013],
[102.9, 0.146, 0.016, 0.074, 0.012],
[119.5, 0.172, 0.017, 0.052, 0.010]]
)

LO_SRP100_475 = np.array(
[[102.9, 0.132, 0.011, 0.077, 0.008],
[119.5, 0.194, 0.013, 0.050, 0.007]]
)

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

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 = np.array(NTPD_475_S + Hdes_475_S)
E_1100 = Hdes_1100_E + Ndes_1100_E
S_1100 = np.array(Hdes_1100_S + Ndes_1100_S)

x_axis = np.linspace(0.,150.,150)
#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.plot(x_axis, np.poly1d(np.polyfit(E_475, S_475, 2))(x_axis), color='b', label='Exp. 475 K')

#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(np.unique(E_1100), np.poly1d(np.polyfit(E_1100, S_1100, 2))(np.unique(E_1100)))
#plt.plot(x_axis, np.poly1d(np.polyfit(E_1100, S_1100, 2))(x_axis), color='r', label='Exp. 1100 K')

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)*1.5 )
#plt.plot(x_axis, fit_475_S, color='b', label=r'Exp. 475 K $\times$ 1.5', 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) )*1.5 )
#plt.plot(x_axis, fit_1100_S, color='r', label=r'Exp. 1100 K $\times$ 1.5', linestyle='--')

#plt.errorbar( LO_SRP32[:,0], LO_SRP32[:,1], xerr=len(LO_SRP32)*[4.2], yerr=LO_SRP32[:,2], marker='o', label='AIMD SRP032', color='k' )
plt.errorbar( LO_SRP100_475[:,0], LO_SRP100_475[:,1], yerr=LO_SRP100_475[:,2], marker='o', label='AIMD 475 K reaction', color='b' )
plt.errorbar( LO_SRP100_475[:,0], LO_SRP100_475[:,3], yerr=LO_SRP100_475[:,4], marker='o', color='b', linestyle='--', fillstyle='none', label='AIMD 475 K trapping' )
plt.errorbar( LO_SRP100[:,0], LO_SRP100[:,1], yerr=LO_SRP100[:,2], marker='o', label='AIMD 1100 K reaction', color='r' )
plt.errorbar( LO_SRP100[:,0], LO_SRP100[:,3], yerr=LO_SRP100[:,4], marker='o', color='r', linestyle='--', fillstyle='none', label='AIMD 1100 K trapping' )
plt.errorbar( LO_1100_FS[:,0], LO_1100_FS[:,1], yerr=LO_1100_FS[:,2], marker='o', label='AIMD 1100 K (FS) reaction', color='k' )
plt.errorbar( LO_1100_FS[:,0], LO_1100_FS[:,3], yerr=LO_1100_FS[:,4], marker='o', color='k', linestyle='--', fillstyle='none', label='AIMD 1100 K (FS) trapping' )

#plt.plot( Guo[:,0], Guo[:,1], marker='o', color='g', label='Guo et al' )
plt.legend(loc='upper left', numpoints=1, frameon=False, fontsize=10)
plt.xlim(60.,130.)
plt.ylim(0,0.3)
plt.tick_params(length=6, width=1)
plt.ylabel('Reaction probability')
plt.xlabel('Incidence energy (kJ/mol)')

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