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

# 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(
[[71.4, 0.022, 0.007],
[89.2, 0.060, 0.011],
[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]]
)

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

plt.subplot(2,1,1)

plt.errorbar( Ni_LO[:,0], Ni_LO[:,1], yerr=Ni_LO[:,2], marker='o', label='Ni(111)', color='b' )
plt.errorbar( LO[:,0], LO[:,1], yerr=LO[:,2], marker='o', label='Pd(111)', color='k' )
plt.errorbar( Pt_LO[:,0], Pt_LO[:,1], yerr=Pt_LO[:,2], marker='o', label='Pt(111)', color='r' )

plt.xlim(50,150.)
plt.ylim(0,0.25)
plt.tick_params(labelbottom='off', length=6, width=1)
plt.ylabel('Reaction probability')
plt.legend(loc='best', numpoints=1)
plt.annotate('(a)', xy=(120, 0.22))

plt.subplot(2,1,2)

plt.errorbar( Ni_v1[:,0], Ni_v1[:,1], yerr=Ni_v1[:,2], marker='o', label='Ni(111)\cite{nattino_chemically_2016}', color='b' )
plt.errorbar( v1[:,0], v1[:,1], yerr=v1[:,2], marker='o', label='Pd(111)', color='k' )
plt.errorbar( Pt_v1[:,0], Pt_v1[:,1], yerr=Pt_v1[:,2], marker='o', label='Pt(111)\cite{migliorini_surface_2017}', color='r' )

plt.xlim(50,150.)
plt.ylim(0,0.28)
plt.tick_params(length=6, width=1)
plt.ylabel('Reaction probability')
plt.xlabel('Incidence energy (kJ/mol)')
plt.annotate('(b)', xy=(120, 0.25))

plt.tight_layout()
plt.subplots_adjust(wspace=0, hspace=0)
plt.savefig('reactionprobability_comparison.pdf')
#plt.show()