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

mpl.rc("text", usetex=True)
fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(3.69,5.5))
plt.rcParams.update({'font.size': 10})

plt.subplot(2,1,1)

plt.errorbar( Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,1], yerr=Shirhatti_DS1[:,2], marker='s', label='Exp. old', color='g', markerfacecolor='None', linestyle='--', capsize=4. )
plt.plot( Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,4], marker='v', label='Exp. new (upper limit)', color='g' )
plt.plot( Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,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.plot( RPBE[:,0]*96.48530749926, RPBE[:,1], marker='d', label='RPBE', color='orange' )
#plt.plot( RPBE_QD[:,0]*96.48530749926, RPBE_QD[:,1], linestyle='--', label='', color='orange' )
plt.errorbar( NN_DS1[:,0]*96.48530749926, NN_DS1[:,1], yerr=NN_DS1[:,2], marker='o', label='MS-RPBEl', color='r', capsize=4., fillstyle='full' )
plt.fill_between(Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,3], Shirhatti_DS1[:,4], color='grey', alpha=0.5)

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

handles, labels = plt.gca().get_legend_handles_labels()
order = [0,1,3,2,4,5]
plt.legend([handles[idx] for idx in order],[labels[idx] for idx in order], loc='best', numpoints=1, frameon=True)

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(60., 270.)
plt.ylim(0.0, 0.6)
#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.subplot(2,1,2)

plt.errorbar( Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,1], yerr=Shirhatti_DS1[:,2], marker='s', label='Exp. old', color='g', markerfacecolor='None', linestyle='--', capsize=4. )
plt.plot( Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,4], marker='v', label='Exp. new (upper limit)', color='g' )
plt.plot( Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,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.plot( RPBE[:,0]*96.48530749926, RPBE[:,1], marker='d', label='RPBE', color='orange' )
plt.plot( RPBE_QD[:,0]*96.48530749926, RPBE_QD[:,1], linestyle='--', label='', color='orange' )
plt.errorbar( NN_DS1[:,0]*96.48530749926, NN_DS1[:,1], yerr=NN_DS1[:,2], marker='o', label='MS-RPBEl', color='r', capsize=4., fillstyle='full' )
plt.fill_between(Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,3], Shirhatti_DS1[:,4], color='grey', alpha=0.5)

dx = 50.
i = 0
plt.annotate('{:4.0f}'.format(dx), xy=(NN_DS1[i,0]*96.48530749926 + 0.2 * dx, NN_DS1[i,1]+0.0003), size=10, color='b')
plt.plot( [NN_DS1[i,0]*96.48530749926 + dx, NN_DS1[i,0]*96.48530749926], [NN_DS1[i,1],NN_DS1[i,1]], zorder=10, c='b' )

dx = 35.
i = 0
plt.annotate('{:4.0f}'.format(dx), xy=(NN_DS1[i,0]*96.48530749926 + 0.2 * dx, NN_DS1[i,1]-0.0006), size=10)
plt.plot( [NN_DS1[i,0]*96.48530749926 + dx, NN_DS1[i,0]*96.48530749926], [NN_DS1[i,1],NN_DS1[i,1]], zorder=10, c='k' )

dx = 107.
i = 2
plt.annotate('{:4.0f}'.format(dx), xy=(NN_DS1[i,0]*96.48530749926 + 0.2 * dx, NN_DS1[i,1]+0.005), size=10, color='b')
plt.plot( [NN_DS1[i,0]*96.48530749926 + dx, NN_DS1[i,0]*96.48530749926], [NN_DS1[i,1],NN_DS1[i,1]], zorder=10, c='b' )

dx = 48.
i = 2
plt.annotate('{:4.0f}'.format(dx), xy=(NN_DS1[i,0]*96.48530749926 + 0.2 * dx, NN_DS1[i,1]-0.017), size=10)
plt.plot( [NN_DS1[i,0]*96.48530749926 + dx, NN_DS1[i,0]*96.48530749926], [NN_DS1[i,1],NN_DS1[i,1]], zorder=10, c='k' )

dx = 72.
i = 4
plt.annotate('{:4.0f}'.format(dx), xy=(NN_DS1[i,0]*96.48530749926 + 0.2 * dx, NN_DS1[i,1]-0.09), size=10)
plt.plot( [NN_DS1[i,0]*96.48530749926 + dx, NN_DS1[i,0]*96.48530749926], [NN_DS1[i,1],NN_DS1[i,1]], zorder=10, c='k' )

i = 6
dy = NN_DS1[i,1] - Shirhatti_DS1[i,3]
factor = NN_DS1[i,1] / Shirhatti_DS1[i,3]
plt.annotate('{:2.1f}x'.format(factor), xy=(NN_DS1[i,0]*96.48530749926 + 5., NN_DS1[i,1]-dy), size=10, color='b')
plt.plot( [NN_DS1[i,0]*96.48530749926, NN_DS1[i,0]*96.48530749926], [NN_DS1[i,1],Shirhatti_DS1[i,3]], zorder=10, c='b' )

i = 6
dy = NN_DS1[i,1] - Shirhatti_DS1[i,4]
factor = NN_DS1[i,1] / Shirhatti_DS1[i,4]
plt.annotate('{:2.1f}x'.format(factor), xy=(NN_DS1[i,0]*96.48530749926 + 5., NN_DS1[i,1]-dy), size=10)
plt.plot( [NN_DS1[i,0]*96.48530749926, NN_DS1[i,0]*96.48530749926], [NN_DS1[i,1],Shirhatti_DS1[i,4]], zorder=10, c='k' )

i = 0
dy = NN_DS1[i,1] - Shirhatti_DS1[i,3]
factor = NN_DS1[i,1] / Shirhatti_DS1[i,3]
plt.annotate('{:2.1f}x'.format(factor), xy=(NN_DS1[i,0]*96.48530749926 - 25., NN_DS1[i,1]-dy), size=10, color='b')
plt.plot( [NN_DS1[i,0]*96.48530749926, NN_DS1[i,0]*96.48530749926], [NN_DS1[i,1],Shirhatti_DS1[i,3]], zorder=10, c='b' )

i = 0
dy = NN_DS1[i,1] - Shirhatti_DS1[i,4]
factor = NN_DS1[i,1] / Shirhatti_DS1[i,4]
plt.annotate('{:2.1f}x'.format(factor), xy=(NN_DS1[i,0]*96.48530749926 - 25., NN_DS1[i,1]-dy + 0.0002), size=10)
plt.plot( [NN_DS1[i,0]*96.48530749926, NN_DS1[i,0]*96.48530749926], [NN_DS1[i,1],Shirhatti_DS1[i,4]], zorder=10, c='k' )

plt.annotate('(b)', xy=(240, 1e-5), size=12)

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(60., 270.)
#plt.ylim(0.0, 0.8)
plt.ylim(1e-6, 0.9)
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.pdf')
#plt.show()