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

mpl.rc("text", usetex=True)
Nrows = 2
Ncols = 3
#fig, ax = plt.subplots(nrows=Nrows, ncols=Ncols, figsize=(6.69,5.))
fig = plt.figure(figsize=(6.69,4.5))
gs1 = fig.add_gridspec(nrows=2, ncols=3, left=0.1, right=0.67, top=0.92, wspace=0., hspace=0.)
gs2 = fig.add_gridspec(nrows=2, ncols=1, left=0.71, right=0.9, top=0.92, wspace=0., hspace=0.)
#plt.rcParams.update({'font.size': 10})

minx = 0.
maxx = 270.

miny0 = -2
maxy0 = 27.5
miny1 = 0.
maxy1 = 0.7

def EnergySticking( E, S, minS=None, maxS=None, Nsteps=None, targetS=None ):
	if targetS == None:
		yaxis = np.linspace(minS, maxS, Nsteps)
	else:
		yaxis = [ targetS ]

	E_S = []
	for i in yaxis:	# Here we determine the energy that S0=i
		yreduced = S - i
		freduced = interpolate.UnivariateSpline(E, yreduced, s=0)
		E_S.append( freduced.roots()[0] )
	E_S = np.array( E_S )
	
	return E_S

#E0 = EnergySticking( NN_DS2_v0_nopara[:,0]*96.48530749926, NN_DS2_v0_nopara[:,1], minS=0.15, maxS=0.3, Nsteps=4 )
#E1 = EnergySticking( NN_DS2_v2_nopara[:,0]*96.48530749926, NN_DS2_v2_nopara[:,1], minS=0.15, maxS=0.3, Nsteps=4 )
#chi = (E0-E1) / E_J[7][2]
#chi = (E0-E1) / E_v[1][1]
#for i in range(len(chi)):
#	print('{0:3.1f}'.format( chi[i] ))
#exit()

#S = min( max( NN_DS2_v0_j0_nopara[:,1] ), max( NN_DS2_v1_j0_nopara[:,1] ) ) - 0.00001
#E0 = EnergySticking( NN_DS2_v0_j0_nopara[:,0]*96.48530749926, NN_DS2_v0_j0_nopara[:,1], targetS=S )
#E1 = EnergySticking( NN_DS2_v1_j0_nopara[:,0]*96.48530749926, NN_DS2_v1_j0_nopara[:,1], targetS=S )
#chi = (E0-E1) / E_J[7][2]
#chi = (E0-E1) / E_v[0][0]
#print('S={:4.3f}, chi={:3.1f}'.format( S, chi[0] ))
#exit()

prop_cycle = plt.rcParams['axes.prop_cycle']
color_cycle = prop_cycle.by_key()['color']

ax1 = fig.add_subplot(gs1[1,0])

minS = max( min(NN_DS1_v0_j2[:,1]), min(NN_DS1_v0_j0[:,1]) ) + 0.00001
maxS = min( max(NN_DS1_v0_j2[:,1]), max(NN_DS1_v0_j0[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS1_v0_j0[:,0]*96.48530749926, NN_DS1_v0_j0[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS1_v0_j2[:,0]*96.48530749926, NN_DS1_v0_j2[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_J[1][0]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax1.plot( xaxis, Efficacy, label=r'$J=0 \rightarrow J=2$', color=color_cycle[1] )

minS = max( min(NN_DS1_v0_j4[:,1]), min(NN_DS1_v0_j0[:,1]) ) + 0.00001
maxS = min( max(NN_DS1_v0_j4[:,1]), max(NN_DS1_v0_j0[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS1_v0_j0[:,0]*96.48530749926, NN_DS1_v0_j0[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS1_v0_j4[:,0]*96.48530749926, NN_DS1_v0_j4[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_J[3][0]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax1.plot( xaxis, Efficacy, label=r'$J=0 \rightarrow J=4$', color=color_cycle[2] )

minS = max( min(NN_DS1_v0_j6[:,1]), min(NN_DS1_v0_j0[:,1]) ) + 0.00001
maxS = min( max(NN_DS1_v0_j6[:,1]), max(NN_DS1_v0_j0[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS1_v0_j0[:,0]*96.48530749926, NN_DS1_v0_j0[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS1_v0_j6[:,0]*96.48530749926, NN_DS1_v0_j6[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_J[5][0]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax1.plot( xaxis, Efficacy, label=r'$J=0 \rightarrow J=6$', color=color_cycle[3] )

minS = max( min(NN_DS1_v0_j8[:,1]), min(NN_DS1_v0_j0[:,1]) ) + 0.00001
maxS = min( max(NN_DS1_v0_j8[:,1]), max(NN_DS1_v0_j0[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS1_v0_j0[:,0]*96.48530749926, NN_DS1_v0_j0[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS1_v0_j8[:,0]*96.48530749926, NN_DS1_v0_j8[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_J[7][0]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax1.plot( xaxis, Efficacy, label=r'$J=0 \rightarrow J=8$', color=color_cycle[4] )

ax1.plot([0.,300.], [1.,1.], color='k', label='', linestyle='--')

ax1.annotate('(e)', xy=(20, 25), size=12)
#plt.annotate(r'$\nu=0$', xy=(30, 22), size=12)

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(minx, maxx)
plt.ylim(miny0, maxy0)
#plt.ylim(1e-4, 0.9)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=True)
plt.ylabel('Rotational efficacy')

ax2 = fig.add_subplot(gs1[0,0])

Order = 3

f = interpolate.UnivariateSpline(NN_DS1_v0_j0[:,0]*96.48530749926, NN_DS1_v0_j0[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS1_v0_j0[:,0]*96.48530749926), max(NN_DS1_v0_j0[:,0]*96.48530749926), 200 )
ax2.plot( xaxis, f(xaxis), label=r'$J=0$' )

f = interpolate.UnivariateSpline(NN_DS1_v0_j2[:,0]*96.48530749926, NN_DS1_v0_j2[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS1_v0_j2[:,0]*96.48530749926), max(NN_DS1_v0_j2[:,0]*96.48530749926), 200 )
ax2.plot( xaxis, f(xaxis), label=r'$J=2$' )

f = interpolate.UnivariateSpline(NN_DS1_v0_j4[:,0]*96.48530749926, NN_DS1_v0_j4[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS1_v0_j4[:,0]*96.48530749926), max(NN_DS1_v0_j4[:,0]*96.48530749926), 200 )
ax2.plot( xaxis, f(xaxis), label=r'$J=4$' )

f = interpolate.UnivariateSpline(NN_DS1_v0_j6[:,0]*96.48530749926, NN_DS1_v0_j6[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS1_v0_j6[:,0]*96.48530749926), max(NN_DS1_v0_j6[:,0]*96.48530749926), 200 )
ax2.plot( xaxis, f(xaxis), label=r'$J=6$' )

f = interpolate.UnivariateSpline(NN_DS1_v0_j8[:,0]*96.48530749926, NN_DS1_v0_j8[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS1_v0_j8[:,0]*96.48530749926), max(NN_DS1_v0_j8[:,0]*96.48530749926), 200 )
ax2.plot( xaxis, f(xaxis), label=r'$J=8$' )

plt.annotate('(a)', xy=(220, 0.03), size=12)

ax2.set_title(r'$\nu=0$', loc='center')

plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(minx, maxx)
plt.ylim(miny1, maxy1)
#plt.ylim(1e-4, 0.55)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False)
plt.ylabel('Reaction probability')
#plt.xlabel('Incidence energy (kJ/mol)')

ax3 = fig.add_subplot(gs1[1,1])

minS = max( min(NN_DS2_v1_j2_nopara[:,1]), min(NN_DS2_v1_j0_nopara[:,1]) ) + 0.00001
maxS = min( max(NN_DS2_v1_j2_nopara[:,1]), max(NN_DS2_v1_j0_nopara[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS2_v1_j0_nopara[:,0]*96.48530749926, NN_DS2_v1_j0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS2_v1_j2_nopara[:,0]*96.48530749926, NN_DS2_v1_j2_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_J[1][1]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax3.plot( xaxis, Efficacy, label=r'$J=0 \rightarrow J=2$', color=color_cycle[1] )

minS = max( min(NN_DS2_v1_j4_nopara[:,1]), min(NN_DS2_v1_j0_nopara[:,1]) ) + 0.00001
maxS = min( max(NN_DS2_v1_j4_nopara[:,1]), max(NN_DS2_v1_j0_nopara[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS2_v1_j0_nopara[:,0]*96.48530749926, NN_DS2_v1_j0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS2_v1_j4_nopara[:,0]*96.48530749926, NN_DS2_v1_j4_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_J[3][1]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax3.plot( xaxis, Efficacy, label=r'$J=0 \rightarrow J=4$', color=color_cycle[2] )

minS = max( min(NN_DS2_v1_j6_nopara[:,1]), min(NN_DS2_v1_j0_nopara[:,1]) ) + 0.00001
maxS = min( max(NN_DS2_v1_j6_nopara[:,1]), max(NN_DS2_v1_j0_nopara[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS2_v1_j0_nopara[:,0]*96.48530749926, NN_DS2_v1_j0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS2_v1_j6_nopara[:,0]*96.48530749926, NN_DS2_v1_j6_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_J[5][1]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax3.plot( xaxis, Efficacy, label=r'$J=0 \rightarrow J=6$', color=color_cycle[3] )

minS = max( min(NN_DS2_v1_j8_nopara[:,1]), min(NN_DS2_v1_j0_nopara[:,1]) ) + 0.00001
maxS = min( max(NN_DS2_v1_j8_nopara[:,1]), max(NN_DS2_v1_j0_nopara[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS2_v1_j0_nopara[:,0]*96.48530749926, NN_DS2_v1_j0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS2_v1_j8_nopara[:,0]*96.48530749926, NN_DS2_v1_j8_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_J[7][1]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax3.plot( xaxis, Efficacy, label=r'$J=0 \rightarrow J=8$', color=color_cycle[4] )

ax3.plot([0.,300.], [1.,1.], color='k', label='', linestyle='--')

ax3.annotate('(f)', xy=(20, 25), size=12)
#plt.annotate(r'$\nu=1$', xy=(30, 22), size=12)

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

ax4 = fig.add_subplot(gs1[0,1])

Order = 3

f = interpolate.UnivariateSpline(NN_DS2_v1_j0_nopara[:,0]*96.48530749926, NN_DS2_v1_j0_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v1_j0_nopara[:,0]*96.48530749926), max(NN_DS2_v1_j0_nopara[:,0]*96.48530749926), 200 )
ax4.plot( xaxis, f(xaxis), label=r'$J=0$' )

f = interpolate.UnivariateSpline(NN_DS2_v1_j2_nopara[:,0]*96.48530749926, NN_DS2_v1_j2_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v1_j2_nopara[:,0]*96.48530749926), max(NN_DS2_v1_j2_nopara[:,0]*96.48530749926), 200 )
ax4.plot( xaxis, f(xaxis), label=r'$J=2$' )

f = interpolate.UnivariateSpline(NN_DS2_v1_j4_nopara[:,0]*96.48530749926, NN_DS2_v1_j4_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v1_j4_nopara[:,0]*96.48530749926), max(NN_DS2_v1_j4_nopara[:,0]*96.48530749926), 200 )
ax4.plot( xaxis, f(xaxis), label=r'$J=4$' )

f = interpolate.UnivariateSpline(NN_DS2_v1_j6_nopara[:,0]*96.48530749926, NN_DS2_v1_j6_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v1_j6_nopara[:,0]*96.48530749926), max(NN_DS2_v1_j6_nopara[:,0]*96.48530749926), 200 )
ax4.plot( xaxis, f(xaxis), label=r'$J=6$' )

f = interpolate.UnivariateSpline(NN_DS2_v1_j8_nopara[:,0]*96.48530749926, NN_DS2_v1_j8_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v1_j8_nopara[:,0]*96.48530749926), max(NN_DS2_v1_j8_nopara[:,0]*96.48530749926), 200 )
ax4.plot( xaxis, f(xaxis), label=r'$J=8$' )

ax4.annotate('(b)', xy=(220, 0.03), size=12)

ax4.set_title(r'$\nu=1$', loc='center')

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(minx, maxx)
plt.ylim(miny1, maxy1)
#plt.ylim(1e-4, 0.55)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelleft=False, labelbottom=False)
#plt.xlabel('Incidence energy (kJ/mol)')

ax5 = fig.add_subplot(gs1[1,2])

#S0 = np.unique( np.concatenate((NN_DS2_v2_j0_nopara[:,1], NN_DS1_v2_j0[:,1])), axis=0 )
#E0 = np.unique( np.concatenate((NN_DS2_v2_j0_nopara[:,0], NN_DS1_v2_j0[:,0])), axis=0 ) * 96.48530749926

#S1 = np.unique( np.concatenate((NN_DS2_v2_j2_nopara[:,1], NN_DS1_v2_j2[:,1])), axis=0 )
#E1 = np.unique( np.concatenate((NN_DS2_v2_j2_nopara[:,0], NN_DS1_v2_j2[:,0])), axis=0 ) * 96.48530749926
#minS = max( min( S1 ), min( S0 ) ) + 0.00001
#maxS = min( max( S1 ), max( S0 ) ) - 0.00001
minS = max( min(NN_DS2_v2_j2_nopara[:,1]), min(NN_DS2_v2_j0_nopara[:,1]) ) + 0.00001
maxS = min( max(NN_DS2_v2_j2_nopara[:,1]), max(NN_DS2_v2_j0_nopara[:,1]) ) - 0.00001
Nsteps = 30
#E0 = EnergySticking( E0, S0, minS=minS, maxS=maxS, Nsteps=Nsteps )
#E1 = EnergySticking( E1, S1, minS=minS, maxS=maxS, Nsteps=Nsteps )
E0 = EnergySticking( NN_DS2_v2_j0_nopara[:,0]*96.48530749926, NN_DS2_v2_j0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS2_v2_j2_nopara[:,0]*96.48530749926, NN_DS2_v2_j2_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_J[1][2]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax5.plot( xaxis, Efficacy, label=r'$J=0 \rightarrow J=2$', color=color_cycle[1] )

minS = max( min(NN_DS2_v2_j4_nopara[:,1]), min(NN_DS2_v2_j0_nopara[:,1]) ) + 0.00001
maxS = min( max(NN_DS2_v2_j4_nopara[:,1]), max(NN_DS2_v2_j0_nopara[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS2_v2_j0_nopara[:,0]*96.48530749926, NN_DS2_v2_j0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS2_v2_j4_nopara[:,0]*96.48530749926, NN_DS2_v2_j4_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_J[3][2]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax5.plot( xaxis, Efficacy, label=r'$J=0 \rightarrow J=4$', color=color_cycle[2] )

minS = max( min(NN_DS2_v2_j6_nopara[:,1]), min(NN_DS2_v2_j0_nopara[:,1]) ) + 0.00001
maxS = min( max(NN_DS2_v2_j6_nopara[:,1]), max(NN_DS2_v2_j0_nopara[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS2_v2_j0_nopara[:,0]*96.48530749926, NN_DS2_v2_j0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS2_v2_j6_nopara[:,0]*96.48530749926, NN_DS2_v2_j6_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_J[5][2]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax5.plot( xaxis, Efficacy, label=r'$J=0 \rightarrow J=6$', color=color_cycle[3] )

minS = max( min(NN_DS2_v2_j8_nopara[:,1]), min(NN_DS2_v2_j0_nopara[:,1]) ) + 0.00001
maxS = min( max(NN_DS2_v2_j8_nopara[:,1]), max(NN_DS2_v2_j0_nopara[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS2_v2_j0_nopara[:,0]*96.48530749926, NN_DS2_v2_j0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS2_v2_j8_nopara[:,0]*96.48530749926, NN_DS2_v2_j8_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_J[7][2]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax5.plot( xaxis, Efficacy, label=r'$J=0 \rightarrow J=8$', color=color_cycle[4] )

ax5.plot([0.,300.], [1.,1.], color='k', label='', linestyle='--')

ax5.annotate('(g)', xy=(20, 25), size=12)
#plt.annotate(r'$\nu=2$', xy=(30, 22), size=12)

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(minx, maxx)
plt.ylim(miny0, maxy0)
#plt.ylim(1e-4, 0.9)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=True, labelleft=False)
#plt.ylabel('Reaction probability')
plt.xlabel('Incidence energy (kJ/mol)')
ax5.xaxis.set_label_coords(0.2,-0.15)

ax6 = fig.add_subplot(gs1[0,2])

Order = 3

#S = np.unique( np.concatenate((NN_DS2_v2_j0_nopara[:,1], NN_DS1_v2_j0[:,1])), axis=0 )
#E = np.unique( np.concatenate((NN_DS2_v2_j0_nopara[:,0], NN_DS1_v2_j0[:,0])), axis=0 ) * 96.48530749926
#f = interpolate.UnivariateSpline(E, S, s=0, k=Order)
#xaxis = np.linspace( min(E), max(E), 200 )
f = interpolate.UnivariateSpline(NN_DS2_v2_j0_nopara[:,0]*96.48530749926, NN_DS2_v2_j0_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v2_j0_nopara[:,0])*96.48530749926, max(NN_DS2_v2_j0_nopara[:,0])*96.48530749926, 200 )
ax6.plot( xaxis, f(xaxis), label=r'$J=0$' )

#S = np.unique( np.concatenate((NN_DS2_v2_j2_nopara[:,1], NN_DS1_v2_j2[:,1])), axis=0 )
#E = np.unique( np.concatenate((NN_DS2_v2_j2_nopara[:,0], NN_DS1_v2_j2[:,0])), axis=0 ) * 96.48530749926
#f = interpolate.UnivariateSpline(E, S, s=0, k=Order)
#xaxis = np.linspace( min(E), max(E), 200 )
f = interpolate.UnivariateSpline(NN_DS2_v2_j2_nopara[:,0]*96.48530749926, NN_DS2_v2_j2_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v2_j2_nopara[:,0])*96.48530749926, max(NN_DS2_v2_j2_nopara[:,0])*96.48530749926, 200 )
ax6.plot( xaxis, f(xaxis), label=r'$J=2$' )

f = interpolate.UnivariateSpline(NN_DS2_v2_j4_nopara[:,0]*96.48530749926, NN_DS2_v2_j4_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v2_j4_nopara[:,0])*96.48530749926, max(NN_DS2_v2_j4_nopara[:,0])*96.48530749926, 200 )
ax6.plot( xaxis, f(xaxis), label=r'$J=4$' )

f = interpolate.UnivariateSpline(NN_DS2_v2_j6_nopara[:,0]*96.48530749926, NN_DS2_v2_j6_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v2_j6_nopara[:,0])*96.48530749926, max(NN_DS2_v2_j6_nopara[:,0])*96.48530749926, 200 )
ax6.plot( xaxis, f(xaxis), label=r'$J=6$' )

f = interpolate.UnivariateSpline(NN_DS2_v2_j8_nopara[:,0]*96.48530749926, NN_DS2_v2_j8_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v2_j8_nopara[:,0])*96.48530749926, max(NN_DS2_v2_j8_nopara[:,0])*96.48530749926, 200 )
ax6.plot( xaxis, f(xaxis), label=r'$J=8$' )

ax6.annotate('(c)', xy=(220, 0.03), size=12)

ax6.set_title(r'$\nu=2$', loc='center')

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(minx, maxx)
plt.ylim(miny1, maxy1)
#plt.ylim(1e-4, 0.55)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelleft=False, labelbottom=False)

ax7 = fig.add_subplot(gs2[1,0])

minS = max( min(NN_DS2_v1_j0_nopara[:,1]), min(NN_DS2_v0_j0_nopara[:,1]) ) + 0.0001
maxS = min( max(NN_DS2_v1_j0_nopara[:,1]), max(NN_DS2_v0_j0_nopara[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS2_v0_j0_nopara[:,0]*96.48530749926, NN_DS2_v0_j0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS2_v1_j0_nopara[:,0]*96.48530749926, NN_DS2_v1_j0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_v[0][0]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax7.plot( xaxis, Efficacy, label=r'$\nu=0 \rightarrow \nu=1$', color=color_cycle[1] )

minS = max( min(NN_DS2_v1_nopara[:,1]), min(NN_DS2_v0_nopara[:,1]) ) + 0.00001
maxS = min( max(NN_DS2_v1_nopara[:,1]), max(NN_DS2_v0_nopara[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS2_v0_nopara[:,0]*96.48530749926, NN_DS2_v0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS2_v1_nopara[:,0]*96.48530749926, NN_DS2_v1_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_v[1][0]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax7.plot( xaxis, Efficacy, linestyle='--', color=color_cycle[1] )

minS = max( min(NN_DS2_v2_j0_nopara[:,1]), min(NN_DS2_v0_j0_nopara[:,1]) ) + 0.00001
maxS = min( max(NN_DS2_v2_j0_nopara[:,1]), max(NN_DS2_v0_j0_nopara[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS2_v0_j0_nopara[:,0]*96.48530749926, NN_DS2_v0_j0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS2_v2_j0_nopara[:,0]*96.48530749926, NN_DS2_v2_j0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_v[0][1]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax7.plot( xaxis, Efficacy, label=r'$\nu=0 \rightarrow \nu=2$', color=color_cycle[2] )

minS = max( min(NN_DS2_v2_nopara[:,1]), min(NN_DS2_v0_nopara[:,1]) ) + 0.00001
maxS = min( max(NN_DS2_v2_nopara[:,1]), max(NN_DS2_v0_nopara[:,1]) ) - 0.00001
Nsteps = 30
E0 = EnergySticking( NN_DS2_v0_nopara[:,0]*96.48530749926, NN_DS2_v0_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
E1 = EnergySticking( NN_DS2_v2_nopara[:,0]*96.48530749926, NN_DS2_v2_nopara[:,1], minS=minS, maxS=maxS, Nsteps=Nsteps )
Efficacy = (E0-E1) / E_v[1][1]
xaxis = np.linspace( min(E1), max(E1), len(E1) )
ax7.plot( xaxis, Efficacy, linestyle='--', color=color_cycle[2] )

#ax7.plot([0.,300.], [1.,1.], color='k', label='', linestyle='--')

ax7.annotate('(h)', xy=(20, 3.05), size=12)
#plt.annotate(r'$\nu=2$', xy=(30, 22), size=12)

#plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(minx, maxx)
plt.ylim(1., 3.25)
#plt.ylim(1e-4, 0.9)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=True, labelleft=False, labelright=True)
ax7.yaxis.set_label_position("right")
plt.ylabel('Vibrational efficacy')

ax8 = fig.add_subplot(gs2[0,0])

Order = 3

f = interpolate.UnivariateSpline(NN_DS2_v0_j0_nopara[:,0]*96.48530749926, NN_DS2_v0_j0_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v0_j0_nopara[:,0])*96.48530749926, max(NN_DS2_v0_j0_nopara[:,0])*96.48530749926, 200 )
ax8.plot( xaxis, f(xaxis), label=r'$\nu=0$', color=color_cycle[0] )

f = interpolate.UnivariateSpline(NN_DS2_v0_nopara[:,0]*96.48530749926, NN_DS2_v0_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v0_nopara[:,0])*96.48530749926, max(NN_DS2_v0_nopara[:,0])*96.48530749926, 200 )
ax8.plot( xaxis, f(xaxis), linestyle='--', color=color_cycle[0] )

f = interpolate.UnivariateSpline(NN_DS2_v1_j0_nopara[:,0]*96.48530749926, NN_DS2_v1_j0_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v1_j0_nopara[:,0])*96.48530749926, max(NN_DS2_v1_j0_nopara[:,0])*96.48530749926, 200 )
ax8.plot( xaxis, f(xaxis), label=r'$\nu=1$', color=color_cycle[1] )

f = interpolate.UnivariateSpline(NN_DS2_v1_nopara[:,0]*96.48530749926, NN_DS2_v1_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v1_nopara[:,0])*96.48530749926, max(NN_DS2_v1_nopara[:,0])*96.48530749926, 200 )
ax8.plot( xaxis, f(xaxis), linestyle='--', color=color_cycle[1] )

f = interpolate.UnivariateSpline(NN_DS2_v2_j0_nopara[:,0]*96.48530749926, NN_DS2_v2_j0_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v2_j0_nopara[:,0])*96.48530749926, max(NN_DS2_v2_j0_nopara[:,0])*96.48530749926, 200 )
ax8.plot( xaxis, f(xaxis), label=r'$\nu=2$', color=color_cycle[2] )

f = interpolate.UnivariateSpline(NN_DS2_v2_nopara[:,0]*96.48530749926, NN_DS2_v2_nopara[:,1], s=0, k=Order)
xaxis = np.linspace( min(NN_DS2_v2_nopara[:,0])*96.48530749926, max(NN_DS2_v2_nopara[:,0])*96.48530749926, 200 )
ax8.plot( xaxis, f(xaxis), linestyle='--', color=color_cycle[2] )

# xaxis = np.linspace( 100., 250., 150 )
# yaxis = []
# from scipy.special import erf
# for i in range( len(xaxis) ):
	# yaxis.append( 0.5*( 1 + erf( (xaxis[i] - 385.9) / 48.24 ) ) )
# ax8.plot( xaxis, yaxis, linestyle=':', color=color_cycle[0] )
# print(yaxis)
# yaxis = []
# for i in range( len(xaxis) ):
	# yaxis.append( 0.5*( 1 + erf( (xaxis[i] - 164) / 48.24 ) ) )
# ax8.plot( xaxis, yaxis, linestyle=':', color=color_cycle[1] )

ax8.annotate('(d)', xy=(220, 0.03), size=12)

ax8.set_title(r'$T_\textrm{rot}=0,506$ K', loc='center')

plt.legend(loc='best', numpoints=1, frameon=True)
plt.xlim(minx, maxx)
plt.ylim(miny1, maxy1)
#plt.ylim(1e-4, 0.55)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelleft=False, labelright=True, labelbottom=False)
#plt.xlabel('Incidence energy (kJ/mol)')
ax8.yaxis.set_label_position("right")
plt.ylabel('Reaction probability')

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