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

exp_gs = np.array(
[[0.020651995821983293, 0.013771069285436477],
[0.03963227312757267, 0.05820703709013464],
[0.07131300305159005, 0.07718398631904488],
[0.09050832531181463, 0.09833058553669072],
[0.10925307719089816, 0.16827395703197023],
[0.10936571978618437, 0.1560747639625617],
[0.11037182296680093, 0.19711378950171032],
[0.1292906588567797, 0.24820386261699467],
[0.14295625371208542, 0.268219939787413],
[0.14684498330841533, 0.29707052450488436],
[0.16013168943412426, 0.358120251090585],
[0.16280951113113606, 0.368112161304196],
[0.18136993876338914, 0.45801784873123474],
[0.18149282159460967, 0.4447096381100619],
[0.1817385872570505, 0.4180932168677165],
[0.19244987404509795, 0.4580608577221619],
[0.22077692670039095, 0.5402410551539106],
[0.2316827779712045, 0.5591373625248324],
[0.23195926434145053, 0.5291938886271939],
[0.42224076842730485, 0.7717069961291909],
[0.42244557314600595, 0.749526645093903],
[0.5695337620578784, 0.8698757859381078],
[0.7134909988326137, 0.879307043234276],
[0.8560632436971356, 0.8887329244065781],
[0.8572229504167783, 0.9131366866692607],
[0.8576325598541801, 0.8687759845986851],
[0.9542850267270158, 0.901313822270465],
[1.0016512380445257, 0.8715531365842669]]
)

exp_excited = np.array(
[[0.054590697769676844, 0.08820964835029743],
[0.08283582853748928, 0.17926198619616174],
[0.09679071005795992, 0.3179483175292359],
[0.1049777786880206, 0.33128878489360347],
[0.2930857926966641, 0.6591908677575932],
[1.2883548036946784, 0.8715569766727422],
[1.5843949044586, 0.9104140639400331],
[2.0205956745243445, 0.9198706658201402]]
)

gs_carter = np.array(
[[0.05024492740782016, 0.006297229219143663],
[0.10047593129374843, 0.026448362720403296],
[0.14959052187891608, 0.15743073047858935],
[0.1997202637874511, 0.27833753148614626],
[0.24982722175107275, 0.42191435768261976],
[0.29893928078680565, 0.5554156171284635],
[0.3490740857942104, 0.6712846347607054],
[0.39923040897181117, 0.7657430730478592],
[0.4494069845448908, 0.8400503778337532],
[0.4986101793602774, 0.8828715365239295],
[0.5990645924838299, 0.9307304785894207],
[0.6985443590749721, 0.9483627204030227],
[0.8000430363403924, 0.9571788413098237],
[0.8975330050757567, 0.9546599496221664],
[0.9980253914408322, 0.964735516372796]]
)

gs_MS_RPBEl = np.array(
[[0.025, 0.33100],
[0.05, 0.77400],
[0.075, 0.95000],
[0.1, 0.99300]]
)

TN_Trot_MS_RPBEl = np.array(
[[0.01, 0.00100],
[0.025, 0.30133],
[0.05, 0.76467],
[0.075, 0.94333],
[0.1, 0.99100]]
)

gs_HSE03_13X = np.array(
[[0.05, 0.03100],
[0.075, 0.15300],
[0.1, 0.30300],
[0.15, 0.61600],
[0.2, 0.880],
[0.3, 0.99000]]
)

# Nozzle temperature is 300 K and Trot = 0.8*TN
#TN_HSE03_13X = np.array(
#[[0.05, 0.01900],
#[0.075, 0.10400],
#[0.1, 0.24800],
#[0.15, 0.61900],
#[0.2, 0.89800],
#[0.3, 0.99400]]
#)

# Nozzle temperature is 300 K and Trot = 0.05*TN
#TN_Trot_HSE03_13X = np.array(
#[[0.05, 0.02200],
#[0.075, 0.11600],
#[0.1, 0.26800],
#[0.15, 0.62400],
#[0.2, 0.89400],
#[0.3, 0.995]]
#)

# Nozzle temperature is 300 K and Trot = 0.03*TN = 9 K
TN_Trot_HSE03_13X = np.array(
[[0.05, 0.02200],
[0.075, 0.11550],
[0.1, 0.25250],
[0.125, 0.44350],
[0.15, 0.61800],
[0.175, 0.77467],
[0.2, 0.89600],
[0.225, 0.96000],
[0.25, 0.98100]]
)

mpl.rc("text", usetex=True)
fig = plt.figure(figsize=(3.69,3))

#plt.errorbar( Energy_LO, Stick_LO2, yerr=Stick_LO2_error, marker='o', linestyle='--', color='blue', fillstyle='none' )
#plt.errorbar( Energy_LO, Stick_LO, yerr=Stick_LO_error, marker='o', label='Laser off', color='blue' )
#plt.errorbar( Energy_v1, Stick_v12, yerr=Stick_v12_error, marker='o', linestyle='--', color='red', fillstyle='none' )
#plt.errorbar( Energy_v1, Stick_v1, yerr=Stick_v1_error, marker='o', label=r'$\nu_1=1$', color='red' )
#plt.scatter( exp_gs[:,0]*96.4853075, exp_gs[:,1], c='k', marker='d', label='Exp.')
#plt.plot( [0.,150.], [1.,1.], c='g', label='RPBE')
#plt.plot( TN_Trot_MS_RPBEl[:,0]*96.4853075, TN_Trot_MS_RPBEl[:,1], c='b', marker='o', label='MS-RPBEl' )
plt.plot( gs_HSE03_13X[:,0]*96.4853075, gs_HSE03_13X[:,1], c='b', marker='o', label=r'$T_\textrm{rot}=0$ K' )
plt.plot( TN_Trot_HSE03_13X[:,0]*96.4853075, TN_Trot_HSE03_13X[:,1], c='r', marker='o', label=r'$T_\textrm{rot}=9$ K' )
#plt.plot( gs_carter[:,0]*96.4853075, gs_carter[:,1], c='orange', marker='s', label='ECW')
#plt.plot( gs_HSE03_13X[:,0], gs_HSE03_13X[:,1], c='r', marker='o', label='HSE03-1/3X GS' )
#plt.plot( TN_HSE03_13X[:,0], TN_HSE03_13X[:,1], c='orange', marker='o', label='HSE03-1/3X TN, Tr=240' )

plt.legend(loc='lower right', numpoints=1, scatterpoints=1)
plt.xlim(0.,30.)
plt.ylim(0.,1.1)
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.savefig('reactionprobability_rotational.pdf')
#plt.show()