#!/usr/bin/env python
import matplotlib as mpl
mpl.use('pgf')
import matplotlib.pyplot as plt
import numpy as np
from math import erf, exp

# E, reaction probability, error, reaction probability with trapping, error
LO = np.array(
[[89.2, 0.011, 0.003, 0.013, 0.004],
[102.4, 0.022, 0.005, 0.023, 0.005],
[111.8, 0.034, 0.006, 0.034, 0.006],
[120.0, 0.066, 0.011, 0.066, 0.011]]
)

# Trapping probabilities still need to be computed
v1 = np.array(
[[89.2, 0.054, 0.010],
[97.4, 0.100, 0.013],
[102.4, 0.115, 0.010],
[111.8, 0.114, 0.010],
[120.0, 0.162, 0.016]]
)

Ni_LO = np.array(
[[101.1, 0.014, 0.004],
[112.3, 0.028, 0.006],
[121.2, 0.040, 0.007]]
#[130.7, 0.108, 0.010],
#[136.4, 0.124, 0.015],
#[146.6, 0.212, 0.018],
#[160.4, 0.280, 0.021]]
)
Pt_LO = np.array(
[[81.7, 0.023, 0.005],
[89.3, 0.030, 0.005],
[97.4, 0.060, 0.011],
[102.5, 0.086, 0.013],
[111.9, 0.094, 0.013],
[120.1, 0.140, 0.016]]
)
Ni_v1 = np.array(
[[101.1, 0.090, 0.013],
[112.3, 0.104, 0.014],
[121.2, 0.168, 0.017],
[130.7, 0.204, 0.019],
[136.4, 0.248, 0.020],
[146.6, 0.322, 0.021],
[160.4, 0.410, 0.022]]
)
Pt_v1 = np.array(
[[60.7, 0.025, 0.005],
[71.4, 0.048, 0.010],
[81.9, 0.066, 0.011],
[92.2, 0.084, 0.012],
[104.6, 0.138, 0.015]]
)

Pt_oakes = np.array(
[[44.96088235682149, 0.0012768001752754142],
[55.0202778474082, 0.00388256113462756],
[64.95270261355986, 0.009348634542598058],
[74.93457524704147, 0.014816317186064687],
[84.86531129200367, 0.033451978646596096]]
)

Pt_oakes_150 = np.array(
[[55.02391915247305, 0.0004460483135091043],
[64.98273018721059, 0.002935651693741432],
[74.91051097009117, 0.005922110181689568],
[84.85728986635354, 0.024641504273099923]]
)

Pt_luntz = np.array(
[[30.24561403508772, 0.0002835284024147477],
[33.719298245614034, 0.00045249484640932816],
[36.070175438596486, 0.0006762545437171885],
[40.877192982456144, 0.001258153009093861],
[45.789473684210535, 0.0023579063632698572],
[49.12280701754385, 0.003790745160956509],
[49.12280701754385, 0.003575425752237599],
[51.08771929824562, 0.004354646822208203]]
)

Pt_higgins = np.array(
[[4.456140350877186, 0.000004505740352001791],
[5.368421052631572, 0.0000062145162482109705],
[7.017543859649127, 0.000009220640909972132],
[9.999999999999996, 0.000020004780331863742],
[10.315789473684209, 0.000014932596786930868],
[11.824561403508774, 0.00003410335234435916],
[18.771929824561408, 0.00006924979730991275],
[22.10526315789474, 0.000062950255005677],
[23.578947368421048, 0.00007446402938691667],
[25.263157894736842, 0.00010808736099786844],
[27.964912280701746, 0.00018157213922627058],
[29.157894736842103, 0.00014159848243475075],
[30.210526315789473, 0.0001750146293393456],
[32.491228070175445, 0.00015795248938262752],
[35.96491228070175, 0.00012867081427239365],
[37.47368421052632, 0.00017488193810140772],
[37.614035087719294, 0.0002193568954451898],
[43.50877192982454, 0.00031364582539200093],
[45.15789473684213, 0.0002872562477333421]]
)

Pt_beck = np.array(
[[9.684210526315791, 6.303195926129284e-7],
[22.175438596491226, 0.000003885711043532239],
[32.9122807017544, 0.000092632065457812],
[43.47368421052634, 0.00039055229842073914],
[54.14035087719297, 0.0015759580585228059],
[63.68421052631582, 0.0094382093220848]]
)

Pt_LO_120K_nattino_2014 = np.array(
[[49.4621315514925,	0.000187874,	0.000055],
[69.88430822225,	0.00785793,	0.003085],
[75.15240601175,	0.0144911,	0.002951],
[79.88018607925,	0.0279144,	0.002042],
[84.13518814,	0.0325766,	0.001960]]
)

Pt_LO_exp_jpcl_2017 = np.array(
[[81.7, 0.023, 0.005],
[89.3, 0.036, 0.005],
[97.4, 0.054, 0.005],
[102.5, 0.071, 0.006],
[111.9, 0.100, 0.008],
[120.1, 0.130, 0.010]]
)

Ni_beck = np.array(
[[50.19336171051138, 0.0000021232578532663575],
[56.70460671131082, 0.000005173061718975147],
[62.15426797161555, 0.0000157500845588202],
[72.29695261318791, 0.00011789716828355196],
[76.58648270103626, 0.00025587312728840724],
[81.27950608230383, 0.0006290858715916917],
[91.02100377654403, 0.002457186959678674]]
)

Pd_zero = np.array(
[[60.43555181820967, 0.006317090344236713],
[60.714670319042746, 0.0074241857969656365],
[65.72173883171266, 0.010672455243688116],
[70.96315486427932, 0.022779190710074344],
[71.24234804564773, 0.02707103450191421],
[71.24391633688987, 0.03420078129770938]]
)

Pd_twenty = np.array(
[[47.632208818725694, 0.0008167540589646499],
[47.81820069198984, 0.0008978528727679341],
[51.80770957038528, 0.0017139966307244843],
[55.93578818210401, 0.0030606895909412336],
[56.170845167086505, 0.005051362017160901]]
)

Pd_exp = []
for i in np.arange(0., 100., 0.5):
	Pd_exp.append( [i, 1.07 * 10**(-6) * exp( 0.143 * i ) ] )
Pd_exp = np.array( Pd_exp )

mpl.rc("text")
fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(5,5.5))

plt.subplot(2,1,1)

plt.scatter( Ni_beck[:,0], Ni_beck[:,1], marker='s', label=r'Ni(111) Bisson et al.$\cite{bisson_state-resolved_2007}$', color='b' )
plt.plot( Pd_exp[:,0], Pd_exp[:,1], color='k' )
plt.scatter( Pd_zero[:,0], Pd_zero[:,1], marker='o', label='Pd(111) Tait et al. (0$^\circ$)$\cite{tait_methane_2005}$', color='k' )
plt.scatter( Pd_twenty[:,0], Pd_twenty[:,1], marker='v', label='Pd(111) Tait et al. (28$^\circ$)$\cite{tait_methane_2005}$', color='k' )
plt.scatter( Pt_luntz[:,0], Pt_luntz[:,1], marker='^', label='Pt(111) Luntz et al.$\cite{luntz_activation_1989}$', color='r' )
plt.scatter( Pt_oakes[:,0], Pt_oakes[:,1], marker='D', label='Pt(111) Oakes et al.$\cite{oakes_dissociative_1993}$', color='r' )
#plt.scatter( Pt_higgins[:,0], Pt_higgins[:,1], marker='o', label='Pt(111) Higgins et al.$\cite{higgins_state_2001}$', color='r' )
plt.scatter( Pt_beck[:,0], Pt_beck[:,1], marker='v', label='Pt(111) Bisson et al.$\cite{bisson_state-resolved_2007}$', color='r' )

plt.legend(loc='best', scatterpoints=1, frameon=False, fontsize='small', markerfirst=False, handletextpad=-1.5)
plt.xlim(20.,180.)
plt.ylim(10**-6,0.1)
plt.yscale('log')
plt.tick_params(labelbottom='off', length=6, width=1)
plt.ylabel('Reaction probability')
plt.annotate('(a)', xy=(40, 0.01))

plt.subplot(2,1,2)

#plt.scatter( Ni_beck[:,0], Ni_beck[:,1], marker='s', color='b' )
plt.plot( Pd_exp[:,0], Pd_exp[:,1], color='k' )
plt.scatter( Pd_zero[:,0], Pd_zero[:,1], marker='o', label='CH$_4$ + Pd(111) Tait et al. (0$^\circ$)$\cite{tait_methane_2005}$', color='k' )
plt.scatter( Pd_twenty[:,0], Pd_twenty[:,1], marker='v', label='CH$_4$ + Pd(111) Tait et al. (28$^\circ$)$\cite{tait_methane_2005}$', color='k' )
#plt.scatter( Pt_luntz[:,0], Pt_luntz[:,1], marker='^', color='r' )
#plt.scatter( Pt_oakes[:,0], Pt_oakes[:,1], marker='D', color='r' )
#plt.scatter( Pt_higgins[:,0], Pt_higgins[:,1], marker='o', color='r' )

#plt.errorbar( Ni_LO[:,0], Ni_LO[:,1], yerr=Ni_LO[:,2], marker='o', label='Ni(111)', color='b', fillstyle='none' )
plt.errorbar( Pt_LO[:,0], Pt_LO[:,1], yerr=Pt_LO[:,2], marker='o', label='CHD$_3$ + Pt(111) BOMD', color='r', fillstyle='none' )
plt.errorbar( LO[:,0], LO[:,1], yerr=LO[:,2], marker='o', label='CHD$_3$ + Pd(111) BOMD', color='k', fillstyle='none' )

plt.scatter( Pt_oakes[:,0], Pt_oakes[:,1], marker='D', label='CH$_4$ + Pt(111) Oakes et al.$\cite{oakes_dissociative_1993}$', color='r' )
plt.scatter( Pt_LO_120K_nattino_2014[:,0], Pt_LO_120K_nattino_2014[:,1], marker='s', color='r', label='CHD$_3$ + Pt(111) Nattino et al.$\cite{nattino_ab_2014}$' )
plt.scatter( Pt_LO_exp_jpcl_2017[:,0], Pt_LO_exp_jpcl_2017[:,1], marker='d', color='r', label='CHD$_3$ + Pt(111) Migliorini et al.$\cite{migliorini_surface_2017}$' )

plt.legend(loc=(0.34,0.05), scatterpoints=1, numpoints=1, frameon=False, fontsize='small', markerfirst=True)
plt.xlim(20.,180.)
plt.ylim(10**-6,0.4)
plt.yscale('log')
plt.tick_params(length=6, width=1)
plt.ylabel('Reaction probability')
plt.xlabel('Incidence energy (kJ/mol)')
plt.annotate('(b)', xy=(40, 0.05))

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