#!/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, ax = plt.subplots(nrows=2, ncols=1, figsize=(3.69,4.))
plt.rcParams.update({'font.size': 10})

#Ei, Reaction, error
NN_ds = np.array(
[[0.94,0.002491,0.001759],#still running
[1.18,0.031950,0.001244],
[1.29,0.046900,0.001495],
[1.55,0.127350,0.002357],
[1.80,0.195260,0.002803],
[2.12,0.267403,0.003130],
[2.56,0.343756,0.003359]]
)

NN_fs = np.array(
[[0.94, 0.001000, 0.000224],
[1.18, 0.025300, 0.001110],
[1.29, 0.039050, 0.001370],
[1.55, 0.113900, 0.002246],
[1.80, 0.182650, 0.002732],
[2.12, 0.255388, 0.003084],
[2.56, 0.335300, 0.003338]]
)

NN = np.array(
[[0.94, 0.001150, 0.000240],
[1.18, 0.019303, 0.000973],
[1.29, 0.033152, 0.001266],
[1.55, 0.105661, 0.002174],
[1.80, 0.171730, 0.002668],
[2.12, 0.245175, 0.003046],
[2.56, 0.335325, 0.003345]]
)

NN_v2 = np.array(
[[0.41, 0.027999, 0.002123],
[0.61, 0.136348, 0.003704],
[0.82, 0.223456, 0.004148], 
[0.94, 0.260376, 0.003103],
[1.18, 0.313047, 0.003279],
[1.29, 0.339317, 0.003348],
[1.55, 0.392707, 0.003454],
[1.80, 0.429594, 0.003502],
[2.12, 0.466750, 0.003533],
[2.56, 0.511784, 0.003547]]
)

# Ei, S0 and error from 2016 paper, lower and upper limit from Jan Geweke's PhD thesis
Shirhatti = np.array(
[[0.94, 0.000006, 0.000005, 2.4E-5, 6.1E-5],
[1.18, 0.00003, 0.00001, 7.1E-5, 1.8E-4],
[1.29, 0.00012, 0.00007, 3.0E-4, 7.8E-4],
[1.55, 0.0012, 0.0001, 3.1E-3, 8.4E-3],
[1.80, 0.0045, 0.0004, 1.2E-2, 3.2E-2],
[2.12, 0.0082, 0.0008, 2.2E-2, 6.0E-2],
[2.56, 0.021, 0.007, 5.6E-2, 1.6E-1]]
)

SRP32vdW = np.array(
[[1.29, 0.162, 0.017],
[1.80, 0.266, 0.020],
[2.56, 0.382, 0.022]]
)

plt.subplot(2,1,1)

#plt.errorbar( Shirhatti[:,0]*96.48530749926, Shirhatti[:,1], yerr=Shirhatti[:,2], marker='s', label='Exp. old', color='g', markerfacecolor='None', linestyle='--', capsize=4. )
#plt.plot( Shirhatti[:,0]*96.48530749926, Shirhatti[:,4], marker='v', label='Exp. new (upper limit)', color='g' )
#plt.plot( Shirhatti[:,0]*96.48530749926, Shirhatti[:,3], marker='^', label='Exp. new (lower limit)', color='g' )
#plt.errorbar( SRP32vdW[:,0]*96.48530749926, SRP32vdW[:,1], yerr=SRP32vdW[:,2], marker='d', label='SRP32-vdW-DF1', color='k', capsize=4. )
plt.errorbar( NN_fs[:,0]*96.48530749926, NN_fs[:,1], yerr=NN_fs[:,2], marker='o', label='Frozen surface', color='b', capsize=4. )
plt.errorbar( NN_ds[:,0]*96.48530749926, NN_ds[:,1], yerr=NN_ds[:,2], marker='o', label='Distorted surface', color='dimgray', capsize=4. )
plt.errorbar( NN[:,0]*96.48530749926, NN[:,1], yerr=NN[:,2], marker='o', label='Mobile surface', color='r', capsize=4. )
#plt.errorbar( NN_v2[:,0]*96.48530749926, NN_v2[:,1], yerr=NN_v2[:,2], marker='o', label=r'MS.RPBE ($\nu=2,J=1$)', color='orange', capsize=4. )
#plt.fill_between(Shirhatti[:,0]*96.48530749926, Shirhatti[:,3], Shirhatti[:,4], color='grey', alpha=0.5)

plt.annotate('(a)', xy=(240, 0.05), size=12)

plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(70., 270.)
plt.ylim(0.0, 0.4)
#plt.ylim(1e-7, 1.)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False)
plt.ylabel('Reaction probability')
#plt.xlabel('Incidence energy (kJ/mol)')

plt.subplot(2,1,2)

#plt.errorbar( Shirhatti[:,0]*96.48530749926, Shirhatti[:,1], yerr=Shirhatti[:,2], marker='s', label='Exp. old', color='g', markerfacecolor='None', linestyle='--', capsize=4. )
#plt.plot( Shirhatti[:,0]*96.48530749926, Shirhatti[:,4], marker='v', label='Exp. new (upper limit)', color='g' )
#plt.plot( Shirhatti[:,0]*96.48530749926, Shirhatti[:,3], marker='^', label='Exp. new (lower limit)', color='g' )
#plt.errorbar( SRP32vdW[:,0]*96.48530749926, SRP32vdW[:,1], yerr=SRP32vdW[:,2], marker='d', label='SRP32-vdW-DF1', color='k', capsize=4. )
plt.errorbar( NN_fs[:,0]*96.48530749926, NN_fs[:,1], yerr=NN_fs[:,2], marker='o', label='Frozen surface', color='b', capsize=4. )
plt.errorbar( NN_ds[:,0]*96.48530749926, NN_ds[:,1], yerr=NN_ds[:,2], marker='o', label='Distorted surface', color='dimgray', capsize=4. )
plt.errorbar( NN[:,0]*96.48530749926, NN[:,1], yerr=NN[:,2], marker='o', label='Mobile surface', color='r', capsize=4. )
#plt.errorbar( NN_v2[:,0]*96.48530749926, NN_v2[:,1], yerr=NN_v2[:,2], marker='o', label=r'MS.RPBE ($\nu=2,J=1$)', color='orange', capsize=4. )
#plt.fill_between(Shirhatti[:,0]*96.48530749926, Shirhatti[:,3], Shirhatti[:,4], color='grey', alpha=0.5)

plt.annotate('(b)', xy=(240, 2e-3), size=12)

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(70., 270.)
#plt.ylim(0.0, 0.4)
plt.ylim(5e-4, 0.7)
plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True)
plt.ylabel('Reaction probability')
plt.xlabel('Incidence energy (kJ/mol)')

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