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

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

plt.errorbar( NN_DS1[:,0]*96.48530749926, NN_DS1[:,1], yerr=NN_DS1[:,2], marker='d', label=r'DS1 $T_\textrm{vib}=296-1060$ K, $T_\textrm{rot}=11-506$ K', color='C3', capsize=4., fillstyle='full' )
plt.errorbar( NN_DS2[:,0]*96.48530749926, NN_DS2[:,1], yerr=NN_DS2[:,2], marker='s', label=r'DS2 $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K', color='C0', capsize=4., fillstyle='full' )

plt.plot( Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,4], marker='d', label=r'DS1 $T_\textrm{vib}=296-1060$ K, $T_\textrm{rot}=11-506$ K (Exp.)', color='C3', linestyle='--', markerfacecolor='None' )
plt.plot( Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,3], marker='d', label='', color='C3', linestyle='--', markerfacecolor='None' )
plt.fill_between(Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,3], Shirhatti_DS1[:,4], color='C3', alpha=0.3)

plt.plot( Shirhatti_DS2[:,0]*96.48530749926, Shirhatti_DS2[:,4], marker='s', label=r'DS2 $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K (Exp.)', color='C0', linestyle='--', markerfacecolor='None' )
plt.plot( Shirhatti_DS2[:,0]*96.48530749926, Shirhatti_DS2[:,3], marker='s', label='', color='C0', linestyle='--', markerfacecolor='None' )
plt.fill_between(Shirhatti_DS2[:,0]*96.48530749926, Shirhatti_DS2[:,3], Shirhatti_DS2[:,4], color='C0', alpha=0.3)

plt.annotate('(a)', xy=(0.05, 0.8), xycoords='axes fraction', size=12)

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(40., 270.)
plt.ylim(0.0, 0.35)
#plt.ylim(1e-4, 0.9)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False)
#plt.ylabel('Reaction probability')

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

plt.errorbar( NN_DS1[:,0]*96.48530749926, NN_DS1[:,1], yerr=NN_DS1[:,2], marker='d', label=r'Normal $T_\textrm{vib}=296-1060$ K, $T_\textrm{rot}=11-506$ K', color='C3', capsize=4., fillstyle='full' )
plt.errorbar( NN_DS2[:,0]*96.48530749926, NN_DS2[:,1], yerr=NN_DS2[:,2], marker='s', label=r'Off-normal $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K', color='C0', capsize=4., fillstyle='full' )

plt.plot( Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,4], marker='d', label=r'Normal $T_\textrm{vib}=296-1060$ K, $T_\textrm{rot}=11-506$ K (Exp.)', color='C3', linestyle='--', markerfacecolor='None' )
plt.plot( Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,3], marker='d', label='', color='C3', linestyle='--', markerfacecolor='None' )
plt.fill_between(Shirhatti_DS1[:,0]*96.48530749926, Shirhatti_DS1[:,3], Shirhatti_DS1[:,4], color='C3', alpha=0.3)

plt.plot( Shirhatti_DS2[:,0]*96.48530749926, Shirhatti_DS2[:,4], marker='s', label=r'Off-normal $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K (Exp.)', color='C0', linestyle='--', markerfacecolor='None' )
plt.plot( Shirhatti_DS2[:,0]*96.48530749926, Shirhatti_DS2[:,3], marker='s', label='', color='C0', linestyle='--', markerfacecolor='None' )
plt.fill_between(Shirhatti_DS2[:,0]*96.48530749926, Shirhatti_DS2[:,3], Shirhatti_DS2[:,4], color='C0', alpha=0.3)

plt.annotate('(b)', xy=(0.05, 0.8), xycoords='axes fraction', size=12)

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(40., 270.)
#plt.ylim(0.0, 0.55)
plt.ylim(5e-4, 0.55)
plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False, labelleft=False, labelright=True)
plt.tick_params(which='minor', direction='in', axis='y', right=True)
ax2.yaxis.set_label_position("right")
#plt.ylabel('Reaction probability')

ax2 = plt.subplot(Nrows,Ncols,3)

plt.errorbar( NN_DS1[:,0]*96.48530749926, NN_DS1[:,1], yerr=NN_DS1[:,2], marker='d', label=r'DS1 $T_\textrm{vib}=296-1060$ K, $T_\textrm{rot}=11-506$ K', color='C3', capsize=4., fillstyle='full' )
plt.errorbar( NN_DS1_1060[:,0]*96.48530749926, NN_DS1_1060[:,1], yerr=NN_DS1_1060[:,2], marker='d', label=r'DS1 $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K', color='C1', capsize=4., fillstyle='full' )

plt.annotate('(c)', xy=(0.05, 0.8), xycoords='axes fraction', size=12)

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(40., 270.)
plt.ylim(0.0, 0.35)
#plt.ylim(1e-4, 0.9)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False)
plt.ylabel('Reaction probability')
ax2.yaxis.set_label_coords(-0.1,0.)

ax2 = plt.subplot(Nrows,Ncols,4)

plt.errorbar( NN_DS1[:,0]*96.48530749926, NN_DS1[:,1], yerr=NN_DS1[:,2], marker='d', label=r'DS1 $T_\textrm{vib}=296-1060$ K, $T_\textrm{rot}=11-506$ K', color='C3', capsize=4., fillstyle='full' )
plt.errorbar( NN_DS1_1060[:,0]*96.48530749926, NN_DS1_1060[:,1], yerr=NN_DS1_1060[:,2], marker='d', label=r'DS1 $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K', color='C1', capsize=4., fillstyle='full' )

plt.annotate('(d)', xy=(0.05, 0.8), xycoords='axes fraction', size=12)

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(40., 270.)
#plt.ylim(0.0, 0.55)
plt.ylim(5e-4, 0.55)
plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False, labelleft=False, labelright=True)
plt.tick_params(which='minor', direction='in', axis='y', right=True)
ax2.yaxis.set_label_position("right")
ax2.yaxis.set_label_coords(1.14,0.)
plt.ylabel('Reaction probability')

ax2 = plt.subplot(Nrows,Ncols,5)

plt.errorbar( NN_DS1_1060[:,0]*96.48530749926, NN_DS1_1060[:,1], yerr=NN_DS1_1060[:,2], marker='d', label=r'DS1 $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K', color='C1', capsize=4., fillstyle='full' )
plt.errorbar( NN_DS2_nopara[:,0]*96.48530749926, NN_DS2_nopara[:,1], yerr=NN_DS2_nopara[:,2], marker='d', label=r'DS2 (Only normal) $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K', color='C2', capsize=4., fillstyle='full' )

plt.annotate('(e)', xy=(0.05, 0.8), xycoords='axes fraction', size=12)

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

ax2 = plt.subplot(Nrows,Ncols,6)

plt.errorbar( NN_DS1_1060[:,0]*96.48530749926, NN_DS1_1060[:,1], yerr=NN_DS1_1060[:,2], marker='d', label=r'DS1 $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K', color='C1', capsize=4., fillstyle='full' )
plt.errorbar( NN_DS2_nopara[:,0]*96.48530749926, NN_DS2_nopara[:,1], yerr=NN_DS2_nopara[:,2], marker='d', label=r'DS2 (Only normal) $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K', color='C2', capsize=4., fillstyle='full' )

plt.annotate('(f)', xy=(0.05, 0.8), xycoords='axes fraction', size=12)

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(40., 270.)
#plt.ylim(0.0, 0.55)
plt.ylim(5e-4, 0.55)
plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False, labelleft=False, labelright=True)
plt.tick_params(which='minor', direction='in', axis='y', right=True)
ax2.yaxis.set_label_position("right")
#plt.ylabel('Reaction probability')
#plt.xlabel('Normal incidence energy (kJ/mol)')

ax2 = plt.subplot(Nrows,Ncols,7)

plt.errorbar( NN_DS2[:,0]*96.48530749926, NN_DS2[:,1], yerr=NN_DS2[:,2], marker='s', label=r'DS2 $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K', color='C0', capsize=4., fillstyle='full' )
plt.errorbar( NN_DS2_nopara[:,0]*96.48530749926, NN_DS2_nopara[:,1], yerr=NN_DS2_nopara[:,2], marker='d', label=r'DS2 (Only normal) $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K', color='C2', capsize=4., fillstyle='full' )

plt.annotate('(g)', xy=(0.05, 0.8), xycoords='axes fraction', size=12)

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

ax2 = plt.subplot(Nrows,Ncols,8)

plt.errorbar( NN_DS2[:,0]*96.48530749926, NN_DS2[:,1], yerr=NN_DS2[:,2], marker='s', label=r'DS2 $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K', color='C0', capsize=4., fillstyle='full' )
plt.errorbar( NN_DS2_nopara[:,0]*96.48530749926, NN_DS2_nopara[:,1], yerr=NN_DS2_nopara[:,2], marker='d', label=r'DS2 (Only normal) $T_\textrm{vib}=1060$ K, $T_\textrm{rot}=506$ K', color='C2', capsize=4., fillstyle='full' )

plt.annotate('(h)', xy=(0.05, 0.8), xycoords='axes fraction', size=12)

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(40., 270.)
#plt.ylim(0.0, 0.55)
plt.ylim(5e-4, 0.55)
plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelleft=False, labelright=True)
plt.tick_params(which='minor', direction='in', axis='y', right=True)
ax2.yaxis.set_label_position("right")
#plt.ylabel('Reaction probability')
plt.xlabel('Normal incidence energy (kJ/mol)')

lines_labels = [ax.get_legend_handles_labels() for ax in fig.axes]
lines, labels = [sum(lol, []) for lol in zip(*lines_labels)]
lines_ = []
labels_ = []
for i in [0,1,2,3,9,19]:
	lines_.append( lines[i] )
	labels_.append( labels[i] )
#fig.legend([lines[0], lines[1], lines[2]], [labels[0], labels[1], labels[2]], loc='upper center')
fig.legend(lines_, labels_, loc='upper center', ncol=2, fontsize=9, frameon=False)

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