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

#Ei, Reaction, error
NN_111_H2 = np.array(
[[0.025, 0.126500, 0.007433],
[0.035, 0.130000, 0.007520],
[0.043, 0.168000, 0.008360],
[0.132, 0.324000, 0.010465],
[0.181, 0.396000, 0.010936],
[0.233, 0.493000, 0.011179],
[0.282, 0.581500, 0.011031],
[0.338, 0.649000, 0.010672],
[0.413, 0.741871, 0.009788],
[0.454, 0.777500, 0.009300]]
)

NN_211_H2 = np.array(
[[0.025, 0.482500, 0.011161],
[0.035, 0.481000, 0.011167],
[0.043, 0.371000, 0.010644],
[0.132, 0.334500, 0.010633],
[0.181, 0.380000, 0.010963],
[0.233, 0.436500, 0.011150],
[0.282, 0.485000, 0.011180],
[0.338, 0.542271, 0.011074],
[0.413, 0.627884, 0.010570],
[0.454, 0.650445, 0.010908]]
)

Guido_111_H2 = np.array(
[[0.013, 0.113268],
[0.014, 0.111398],
[0.025, 0.120643],
[0.035, 0.122041],
[0.043, 0.149229],
[0.132, 0.269706],
[0.169, 0.300345],
[0.181, 0.338763],
[0.233, 0.431317],
[0.282, 0.525132],
[0.338, 0.621114],
[0.413, 0.71629],
[0.454, 0.754332]]
)

# Guido (111) D2 is not available

Guido_211_H2 = np.array(
[[0.013, 0.707473],
[0.014, 0.679755],
[0.025, 0.586245],
[0.035, 0.588444],
[0.043, 0.464716],
[0.132, 0.336063],
[0.169, 0.34898],
[0.181, 0.346339],
[0.233, 0.378543],
[0.282, 0.422989],
[0.338, 0.482888],
[0.413, 0.566064],
[0.454, 0.61076]]
)

Guido_211_D2 = np.array(
[[0.027, 0.776059],
[0.040, 0.446525],
[0.054, 0.358849],
[0.076, 0.314347],
[0.110, 0.301982],
[0.130, 0.296245],
[0.140, 0.304897],
[0.234, 0.340304],
[0.276, 0.378961],
[0.346, 0.449577],
[0.457, 0.578125]]
)

Groot_211_H2 = np.array(
[[0.00860, 0.51585],
[0.02699, 0.33287],
[0.04029, 0.28664],
[0.05387, 0.27148],
[0.07671, 0.26838],
[0.10991, 0.31772],
[0.13030, 0.31733],
[0.14036, 0.31383],
[0.23418, 0.36084],
[0.27691, 0.42883],
[0.34679, 0.44009],
[0.45694, 0.51158]]
)

AIMD_211_H2 = np.array(
[[0.181, 0.384, 0.022],
[0.338, 0.536, 0.023],
[0.454, 0.651, 0.021]]
)

mpl.rc("text", usetex=True)
Nrows = 3
Ncols = 1
#fig, ax = plt.subplots(nrows=Nrows, ncols=Ncols, figsize=(3.69,4.))
fig = plt.figure(figsize=(3.69, 4.))

outer = mpl.gridspec.GridSpec(2, 1, height_ratios = [1, 2])
gs1 = mpl.gridspec.GridSpecFromSubplotSpec(1, 1, subplot_spec = outer[0])
gs2 = mpl.gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec = outer[1], hspace=0.)

def EnergySticking( E, S, minS=None, maxS=None, Nsteps=None, targetS=None ):
	if targetS == None:
		yaxis = np.linspace(minS, maxS, Nsteps)
	else:
		yaxis = [ targetS ]

	E_S = []
	for i in yaxis:	# Here we determine the energy that S0=i
		yreduced = S - i
		freduced = interpolate.UnivariateSpline(E, yreduced, s=0)
		E_S.append( freduced.roots()[0] )
	E_S = np.array( E_S )
	
	return E_S

ax = plt.subplot(gs1[0])

lines = open('energy_entiredataset.out', 'r').readlines()
Ediff = []
for i in range(1, len(lines) ): Ediff.append( abs( float(lines[i].split()[3]) - float(lines[i].split()[2]) )*2625.5002 )
bins = np.linspace(0.,10.,21)
hist = plt.hist( Ediff, bins=bins, color='b', density=False, weights=np.ones(len(Ediff)) / len(Ediff) )

Tot = 0
for i in range(len(hist[:][0])):
	Tot += hist[0][i]
	print(hist[1][i], Tot)

plt.annotate(r'(a)', xy=(0.12, 0.75), xycoords='axes fraction', size=12)

ylim = plt.ylim()

plt.plot([4.2,4.2], [0.,1.], linestyle='--', color='k')
plt.xlim(0.,10.)
#plt.ylim(0.,0.7)
plt.ylim(0., ylim[1])
plt.xlabel('Absolute error (kJ/mol)')
plt.ylabel('Fraction')

#plt.subplot(Nrows,Ncols,1)
ax = plt.subplot(gs2[0])
plt.plot( Guido_111_H2[:,0]*96.48530749926, Guido_111_H2[:,1], marker='s', label='CRP', mfc='None' )
plt.errorbar( NN_111_H2[:,0]*96.48530749926, NN_111_H2[:,1], yerr=NN_111_H2[:,2], marker='o', label='HD-NNP', capsize=4., mfc='None' )

# for i in range( len( NN_111_H2 ) - 1 ):
	# S = NN_111_H2[i,1]
	# E0 = NN_111_H2[i,0] * 96.48530749926
	# E1 = EnergySticking( Guido_111_H2[:,0], Guido_111_H2[:,1], targetS=S )[0] * 96.48530749926
	# plt.plot( [E0, E1], [S, S], c='k' )
	# print(E1-E0)

plt.annotate(r'(b) H$_2$ + Pt(111)', xy=(0.12, 0.75), xycoords='axes fraction', size=12)

#plt.legend(loc='best', numpoints=1, frameon=True, fontsize=9)
plt.xlim(0., 50.)
plt.ylim(0.0, 0.8)
#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('Sticking probability')
ax.yaxis.set_label_coords(-0.1,0.)

ax = plt.subplot(gs2[1])
#plt.subplot(Nrows,Ncols,2)

plt.plot( Guido_211_H2[:,0]*96.48530749926, Guido_211_H2[:,1], marker='s', label='CRP', linestyle='--', mfc='None' )
plt.errorbar( NN_211_H2[:,0]*96.48530749926, NN_211_H2[:,1], yerr=NN_211_H2[:,2], marker='o', label='HDNNP', capsize=4., mfc='None' )
plt.errorbar( AIMD_211_H2[:,0]*96.48530749926, AIMD_211_H2[:,1], yerr=AIMD_211_H2[:,2], marker='d', label='DFMD', capsize=4., linestyle='', mfc='None' )
#plt.plot( Groot_211_H2[:,0]*96.48530749926, Groot_211_H2[:,1], marker='^', label='Exp.', color='C02', linestyle='--', mfc='None' )

# for i in range( 4, len( NN_211_H2 ) - 2 ):
	# S = NN_211_H2[i,1]
	# E0 = NN_211_H2[i,0] * 96.48530749926
	# E1 = EnergySticking( Guido_211_H2[5:,0], Guido_211_H2[5:,1], targetS=S )[0] * 96.48530749926
	# plt.plot( [E0, E1], [S, S], c='k' )
	# print(E1-E0)

plt.annotate(r'(c) H$_2$ + Pt(211)', xy=(0.12, 0.75), xycoords='axes fraction', size=12)

plt.legend(loc='best', numpoints=1, frameon=True, fontsize=8)
plt.xlim(0., 50.)
plt.ylim(0.0, 0.79)
#plt.ylim(1e-7, 1.)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True)
#plt.ylabel('Sticking probability')
plt.xlabel('Incidence energy (kJ/mol)')

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