#!/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 = plt.figure(figsize=(3.69,3))
#plt.rcParams.update({'font.size': 10})

#Ei, Reaction, error 
Geweke_2016_12 = np.array(
[[0.59, 0.00137049, 0.000662441],
[0.64, 0.00270825, 0.000648565],
[0.92, 0.00738576, 0.00122523],
[0.94, 0.0149115, 0.00244504],
[1.04, 0.0184508, 0.00312071]]
)

Geweke_2016_01 = np.array(
[[0.59, 0.000041764, 0.0000312661],
[0.64, 0.0000924593, 0.0000164302],
[0.86, 0.000639136, 0.000217479],
[0.94, 0.000577473, 0.0000897007],
[1.12, 0.00087208, 0.000276175],
[1.37, 0.00197916, 0.000549222]]
)

#Ei, excitation probability, standard error
# dt = 0.4 fs
NN_12 = np.array(
[[0.28, 0.0020, 0.0002],
[0.32, 0.0031, 0.0003],
[0.45, 0.0087, 0.0004],
[0.52, 0.0133, 0.0005],
[0.78, 0.0287, 0.0008],
[0.94, 0.0342, 0.0009],
[1.27, 0.0275,  0.0009]]
)

# dt = 0.2 fs
NN_12 = np.array(
[[0.28, 0.000194, 0.000044],
[0.32, 0.000419, 0.000065],
[0.45, 0.002074, 0.000145],
[0.52, 0.005293, 0.000232],
[0.78, 0.0230, 0.0005],
[0.97, 0.0295, 0.0006],
[1.27, 0.0284, 0.0006]]
)

# dt = 0.4 fs
NN_12_Geweke = np.array(
[[0.59, 0.013853, 0.000371],
[0.64, 0.016502, 0.000406],
[0.92, 0.031300, 0.000595],
[0.94, 0.032718, 0.000612],
[1.04, 0.034039, 0.000643]]
)

# dt = 0.2 fs
NN_12_Geweke = np.array(
[[0.59, 0.009804, 0.000313],
[0.64, 0.012217, 0.000350],
[0.92, 0.028161, 0.000564],
[0.94, 0.029873, 0.000584],
[1.04, 0.030382, 0.000608]]
)

plt.errorbar( Geweke_2016_12[:,0]*96.48530749926, Geweke_2016_12[:,1], yerr=Geweke_2016_12[:,2], capsize=4, linestyle='', marker='s', label='Exp. (Geweke et al.)', color='k' )
plt.errorbar( NN_12[:,0]*96.48530749926, NN_12[:,1], yerr=NN_12[:,2], marker='o', label='MS-RPBEl (Rahinov et al.)', color='r', capsize=4, linestyle='', fillstyle='full' )
plt.errorbar( NN_12_Geweke[:,0]*96.48530749926, NN_12_Geweke[:,1], yerr=NN_12_Geweke[:,2], marker='o', label='MS-RPBEl (Geweke et al.)', color='b', capsize=4, linestyle='', fillstyle='full' )
#plt.errorbar( NN_12_gauss_nu[:,0]*96.48530749926, NN_12_gauss_nu[:,1], yerr=NN_12_gauss_nu[:,2], marker='o', label=r'MS-RPBEl (GB $\nu$)', color='r', capsize=4, linestyle='', fillstyle='left' )
#plt.errorbar( NN_12_gauss_nuj[:,0]*96.48530749926, NN_12_gauss_nuj[:,1], yerr=NN_12_gauss_nuj[:,2], marker='o', label=r'MS-RPBEl (GB $\nu, j$)', color='r', capsize=4, linestyle='', fillstyle='right' )
#plt.errorbar( NN_12_1GB[:,0]*96.48530749926, NN_12_1GB[:,1], yerr=NN_12_1GB[:,2], marker='o', label=r'MS-RPBEl (1GB $\nu, j, an, am$)', color='r', capsize=4, linestyle='', fillstyle='full' )
#plt.errorbar( NN_12_gauss_nuj[:,0]*96.48530749926, NN_12_gauss_nuj[:,1], yerr=NN_12_gauss_nuj[:,2], marker='o', label=r'MS-RPBEl (GB $\nu, j$)', color='r', capsize=4, linestyle='' )
#plt.errorbar( Gernot_12[:,0]*96.48530749926, Gernot_12[:,1], yerr=Gernot_12[:,2], marker='d', label='SRP32-vdW-DF1 (LDFA)', color='k', capsize=4, linestyle='' )

plt.annotate(r'$\nu=1, j=1\rightarrow\nu=2$', xy=(100., 13*10**-4), size=10)
#plt.annotate('(a)', xy=(240., 5*10**-2), size=12)

plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(0., 200.)
#plt.ylim(0.0, 0.05)
plt.ylim(1e-4, 0.1)
plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True)
plt.ylabel('Vibrational transition probability')
plt.xlabel('Incidence energy (kJ/mol)')

plt.tight_layout()
plt.savefig('Tvjprobability_beamparameters.pdf')
#plt.show()