#!/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)
plt.style.use('tableau-colorblind10')
Nrows = 1
Ncols = 2
fig, ax = plt.subplots(nrows=Nrows, ncols=Ncols, figsize=(6.69,3.))

ax1 = plt.subplot(Nrows,Ncols,1)

plt.errorbar( QCT_LO_CH4_bisson_2007[:,0], QCT_LO_CH4_bisson_2007[:,1], yerr=QCT_LO_CH4_bisson_2007[:,2], marker='o', color='C00', label=r'QCT', capsize=4 )
plt.errorbar( TRPMD_0323K_MS_beta1000K_CH4_bisson_2007[:,0], TRPMD_0323K_MS_beta1000K_CH4_bisson_2007[:,1], yerr=TRPMD_0323K_MS_beta1000K_CH4_bisson_2007[:,2], color='C01', marker='s', label=r'RPMD', capsize=4, ls='--' )
#plt.errorbar( SSRPMD_0323K_MS_beta1000K_CH4_bisson_2007[:,0], SSRPMD_0323K_MS_beta1000K_CH4_bisson_2007[:,1], yerr=SSRPMD_0323K_MS_beta1000K_CH4_bisson_2007[:,2], color='C05', marker='s', label=r'SSRPMD (LO)', capsize=4, ls='--' )
plt.errorbar( Exp_LO_CH4_bisson_2007[:,0], Exp_LO_CH4_bisson_2007[:,1], yerr=Exp_LO_CH4_bisson_2007[:,3], marker='d', color='C02', label=r'Exp.', capsize=4, ls=':' )
plt.plot( RPH_GS_500_CH4[:,1], RPH_GS_500_CH4[:,4], color='C03', marker='^', label=r'RPH', ls='-.' )

#plt.plot([77.2,77.2],[0.,1.], c='k', ls=':')

plt.annotate( r'(A) CH$_4$', xy=(0.1, 0.8), xycoords='axes fraction', fontsize=10)

#plt.title(r'CH$_4$ (laser-off)')
plt.legend(loc='lower right', numpoints=1, fontsize=10)
#plt.legend(loc='best', numpoints=1, borderaxespad=0.4, handletextpad=0.5)
plt.xlim(0,70.)
plt.yscale('log')
plt.ylim(1e-7,10e-2)
locmaj = mpl.ticker.LogLocator(base=10,numticks=12)
ax1.yaxis.set_major_locator(locmaj)
locmin = mpl.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),numticks=12)
ax1.yaxis.set_minor_locator(locmin)
ax1.yaxis.set_minor_formatter(mpl.ticker.NullFormatter())
plt.tick_params(axis='y', which=('minor'), direction='in', right=True)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True)
plt.ylabel('Sticking probability')
plt.xlabel('Incidence energy (kJ/mol)')

ax1 = plt.subplot(Nrows,Ncols,2)

plt.errorbar( QCT_LO_CHD3_migliorini_2017[:,0], QCT_LO_CHD3_migliorini_2017[:,1], yerr=QCT_LO_CHD3_migliorini_2017[:,2], color='C00', marker='o', label=r'QCT', capsize=4 )
plt.errorbar( TRPMD_beta1000K_CHD3_migliorini_2017[:,0], TRPMD_beta1000K_CHD3_migliorini_2017[:,1], yerr=TRPMD_beta1000K_CHD3_migliorini_2017[:,2], color='C01', marker='s', label=r'RPMD', capsize=4, ls='--' )
#plt.errorbar( SSRPMD_beta1000K_CHD3_migliorini_2017[:,0], SSRPMD_beta1000K_CHD3_migliorini_2017[:,1], yerr=SSRPMD_beta1000K_CHD3_migliorini_2017[:,2], color='C05', marker='s', label=r'SSRPMD (LO)', capsize=4, ls='--' )
plt.errorbar( Exp_LO_CHD3_migliorini_2017[:,0], Exp_LO_CHD3_migliorini_2017[:,1], yerr=Exp_LO_CHD3_migliorini_2017[:,2], color='C02', marker='d', label=r'Exp.', capsize=4 )
plt.plot( RPH_LO_500_CHD3[:,0], RPH_LO_500_CHD3[:,1], color='C03', marker='^', label=r'RPH', ls='-.' )

#plt.plot([77.2,77.2],[0.,1.], c='k', ls=':')

plt.annotate( r'(B) CHD$_3$', xy=(0.1, 0.8), xycoords='axes fraction', fontsize=10)

#plt.title(r'CHD$_3$ (laser-off)')
#plt.legend(loc='best', numpoints=1)
#plt.legend(loc='best', numpoints=1, borderaxespad=0.4, handletextpad=0.5)
#plt.xlim(65,125.)
plt.ylim(0,0.2)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelleft=False, labelright=True)
ax1.yaxis.set_label_position("right")
plt.ylabel('Sticking probability')
plt.xlabel('Incidence energy (kJ/mol)')

# ax1 = plt.subplot(Nrows,Ncols,3)

# plt.errorbar( SSRPMD_v3_0_2_MS_beta1000K_CH4_bisson_2007[1:,0], SSRPMD_v3_0_2_MS_beta1000K_CH4_bisson_2007[1:,1] - SSRPMD_0323K_MS_beta1000K_CH4_bisson_2007[:,1] + TRPMD_0323K_MS_beta1000K_CH4_bisson_2007[2:,1], yerr=SSRPMD_v3_0_2_MS_beta1000K_CH4_bisson_2007[1:,2] - SSRPMD_0323K_MS_beta1000K_CH4_bisson_2007[:,2] + TRPMD_0323K_MS_beta1000K_CH4_bisson_2007[2:,2], color='C08', marker='s', label=r'SSRPMDcorr (323 K + $\nu_{3(0)}=2$)', capsize=4, ls='--' )
# plt.errorbar( SSRPMD_v3_0_2_MS_beta1000K_CH4_bisson_2007[:,0], SSRPMD_v3_0_2_MS_beta1000K_CH4_bisson_2007[:,1], yerr=SSRPMD_v3_0_2_MS_beta1000K_CH4_bisson_2007[:,2], color='C05', marker='s', label=r'SSRPMD (323 K + $\nu_{3(0)}=2$)', capsize=4, ls='--' )
# plt.errorbar( SSRPMD_v3_0_1_v3_1_1_MS_beta1000K_CH4_bisson_2007[:,0], SSRPMD_v3_0_1_v3_1_1_MS_beta1000K_CH4_bisson_2007[:,1], yerr=SSRPMD_v3_0_1_v3_1_1_MS_beta1000K_CH4_bisson_2007[:,2], color='C06', marker='s', label=r'SSRPMD (323 K + $\nu_{3(0)}=1 + \nu_{3(1)}=1$)', capsize=4, ls='--' )
# plt.errorbar( SSRPMD_v1_1_v3_0_1_MS_beta1000K_CH4_bisson_2007[:,0], SSRPMD_v1_1_v3_0_1_MS_beta1000K_CH4_bisson_2007[:,1], yerr=SSRPMD_v1_1_v3_0_1_MS_beta1000K_CH4_bisson_2007[:,2], color='C07', marker='s', label=r'SSRPMD (323 K + $\nu_{1}=1 + \nu_{3(0)}=1$)', capsize=4, ls='--' )
# plt.errorbar( QCT_2v3_CH4_bisson_2007[:,0], QCT_2v3_CH4_bisson_2007[:,1], yerr=QCT_2v3_CH4_bisson_2007[:,2], marker='o', color='C00', label=r'QCT ($\nu_{3(0)}=2$)', capsize=4 )
# plt.errorbar( Exp_2v3_CH4_bisson_2007[:,0], Exp_2v3_CH4_bisson_2007[:,1], yerr=Exp_2v3_CH4_bisson_2007[:,3], marker='d', color='C01', label=r'Exp. ($2\nu_{3}\textrm{-F2}$)', capsize=4 )
# plt.plot( RPH_2v3_500K_CH4[:,0], RPH_2v3_500K_CH4[:,1], color='C04', marker='^', label=r'RPH ($2\nu_{3}\textrm{-F2}$)', ls='-.' )

# plt.plot([77.2,77.2],[0.,1.], c='k', ls=':')

# plt.title(r'CH$_4$ (laser-on)')
# plt.legend(loc='best', numpoints=1, fontsize=7)
# plt.xlim(5,70.)
# plt.yscale('log')
# plt.ylim(1e-5,2e-1)
# locmaj = mpl.ticker.LogLocator(base=10,numticks=12)
# ax1.yaxis.set_major_locator(locmaj)
# locmin = mpl.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),numticks=12)
# ax1.yaxis.set_minor_locator(locmin)
# ax1.yaxis.set_minor_formatter(mpl.ticker.NullFormatter())
# plt.tick_params(axis='y', which=('minor'), direction='in', right=True)

# plt.tick_params(length=6, width=1, direction='in', top=True, right=True)
# plt.ylabel('Sticking probability')
# plt.xlabel('Incidence energy (kJ/mol)')

# ax1 = plt.subplot(Nrows,Ncols,4)

# # Temporary RPH results taken from SI in Davide's 2017 JPCL
# tmp = np.array(
# [[60.03569, 0.03287],
# [61.10567, 0.03557],
# [63.29567, 0.04009],
# [64.39072, 0.04226],
# [65.48573, 0.04474],
# [66.58070, 0.04772],
# [67.67575, 0.05028],
# [68.77077, 0.05324],
# [69.86582, 0.05604],
# [70.96089, 0.05891],
# [72.05591, 0.06230],
# [73.15100, 0.06518],
# [74.24603, 0.06894],
# [75.34111, 0.07227],
# [76.43613, 0.07641],
# [77.53120, 0.08028],
# [78.42740, 0.08062],
# [79.77138, 0.08517],
# [81.21492, 0.09025],
# [82.70833, 0.09448],
# [84.15183, 0.10082],
# [85.69511, 0.10427],
# [87.18850, 0.10946],
# [88.68186, 0.11558],
# [89.56361, 0.12005]])
# plt.plot( tmp[:,0], tmp[:,1], color='C04', marker='', label=r'RPH ($\nu_{1}=1$)', ls='-.' )

# plt.errorbar( QCT_v1_CHD3_migliorini_2017[:,0], QCT_v1_CHD3_migliorini_2017[:,1], yerr=QCT_v1_CHD3_migliorini_2017[:,2], color='C00', marker='o', label=r'QCT ($\nu_1=1$)', capsize=4 )
# plt.errorbar( SSRPMD_v1_CHD3_migliorini_2017[:,0], SSRPMD_v1_CHD3_migliorini_2017[:,1], yerr=SSRPMD_v1_CHD3_migliorini_2017[:,2], color='C05', marker='s', label=r'SSRPMD ($\nu_1=1$)', capsize=4 )
# plt.scatter( [SSRPMD_v1_CHD3_migliorini_2017[-1,0]], [SSRPMD_v1_CHD3_migliorini_2017[-1,1] + 0.054415 - 0.009009], label=r'SSRPMDcorr ($\nu_1=1$)', color='C08', marker='s', zorder=10 )
# plt.errorbar( Exp_v1_CHD3_migliorini_2017[:,0], Exp_v1_CHD3_migliorini_2017[:,1], yerr=Exp_v1_CHD3_migliorini_2017[:,2], color='C01', marker='d', label=r'Exp. ($\nu_1=1$)', capsize=4 )

# plt.plot([77.2,77.2],[0.,1.], c='k', ls=':')

# plt.title(r'CHD$_3$ (laser-on)')
# plt.legend(loc='best', numpoints=1)
# #plt.legend(loc='best', numpoints=1, borderaxespad=0.4, handletextpad=0.5)
# #plt.xlim(65,125.)
# plt.ylim(0,0.2)
# plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelleft=False, labelright=True)
# ax1.yaxis.set_label_position("right")
# plt.ylabel('Sticking probability')
# plt.xlabel('Incidence energy (kJ/mol)')

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