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

# J, S0 for theta 10, 30, ..., 170, error
S0_v0 = np.array(
[[0, 0.0000, 0.0019, 0.0132, 0.1494, 0.3149, 0.3544, 0.3224, 0.3913, 0.3926, 0.4464, 0.4814, 0.4561, 0.0000, 0.0019, 0.0040, 0.0112, 0.0131, 0.0132, 0.0129, 0.0141, 0.0159, 0.0182, 0.0249, 0.0466],
[2, 0.3333, 0.3773, 0.2808, 0.3733, 0.2962, 0.2363, 0.2861, 0.2859, 0.2693, 0.2901, 0.3945, 0.3223, 0.0348, 0.0191, 0.0182, 0.0134, 0.0128, 0.0125, 0.0134, 0.0125, 0.0139, 0.0178, 0.0234, 0.0425],
[4, 0.5320, 0.4973, 0.4565, 0.4700, 0.4513, 0.3648, 0.2913, 0.2104, 0.1478, 0.0745, 0.0188, 0.0000, 0.0350, 0.0211, 0.0170, 0.0175, 0.0127, 0.0139, 0.0133, 0.0113, 0.0114, 0.0102, 0.0062, 0.0000],
[6, 0.4468, 0.4037, 0.3911, 0.4099, 0.4374, 0.4344, 0.3788, 0.3240, 0.2309, 0.1517, 0.0954, 0.0455, 0.0363, 0.0202, 0.0174, 0.0150, 0.0137, 0.0142, 0.0138, 0.0136, 0.0136, 0.0126, 0.0163, 0.0199],
[8, 0.2019, 0.2032, 0.2257, 0.2565, 0.2746, 0.3245, 0.3855, 0.4247, 0.5000, 0.5467, 0.5871, 0.6379, 0.0275, 0.0169, 0.0145, 0.0136, 0.0122, 0.0135, 0.0137, 0.0147, 0.0162, 0.0185, 0.0246, 0.0446]]
)

S0_v1 = np.array(
[[0, 0.0000, 0.0000, 0.0327, 0.1857, 0.4057, 0.4968, 0.4798, 0.5634, 0.7256, 0.7260, 0.7364, 0.7748, 0.0000, 0.0000, 0.0062, 0.0123, 0.0141, 0.0140, 0.0141, 0.0147, 0.0149, 0.0169, 0.0224, 0.0396],
[2, 0.5562, 0.4976, 0.4092, 0.5189, 0.4213, 0.3741, 0.4264, 0.4481, 0.4173, 0.4340, 0.6329, 0.7168, 0.0372, 0.0200, 0.0203, 0.0140, 0.0142, 0.0144, 0.0149, 0.0138, 0.0157, 0.0198, 0.0234, 0.0424],
[4, 0.7111, 0.7249, 0.7202, 0.7058, 0.6551, 0.5141, 0.4157, 0.3274, 0.2482, 0.1312, 0.0605, 0.0602, 0.0338, 0.0195, 0.0157, 0.0164, 0.0123, 0.0146, 0.0145, 0.0131, 0.0140, 0.0133, 0.0111, 0.0261],
[6, 0.7088, 0.6837, 0.6974, 0.7261, 0.6604, 0.6288, 0.5178, 0.3798, 0.2778, 0.2101, 0.1957, 0.1518, 0.0337, 0.0195, 0.0166, 0.0138, 0.0133, 0.0141, 0.0144, 0.0143, 0.0146, 0.0143, 0.0221, 0.0339],
[8, 0.5023, 0.5486, 0.5716, 0.5065, 0.4744, 0.5098, 0.5582, 0.5909, 0.6397, 0.6839, 0.6572, 0.6842, 0.0343, 0.0211, 0.0172, 0.0158, 0.0138, 0.0146, 0.0143, 0.0149, 0.0159, 0.0176, 0.0241, 0.0435]]
)

S0_v2 = np.array(
[[0, 0.0048, 0.0141, 0.0503, 0.2008, 0.4619, 0.5459, 0.5474, 0.6561, 0.8207, 0.8576, 0.8716, 0.8416, 0.0048, 0.0053, 0.0077, 0.0128, 0.0142, 0.0141, 0.0144, 0.0147, 0.0134, 0.0137, 0.0175, 0.0363],
[2, 0.5511, 0.5353, 0.4468, 0.5759, 0.4619, 0.4276, 0.4799, 0.4940, 0.4387, 0.5000, 0.6959, 0.8100, 0.0375, 0.0200, 0.0211, 0.0142, 0.0146, 0.0153, 0.0153, 0.0141, 0.0159, 0.0203, 0.0234, 0.0392],
[4, 0.8708, 0.8110, 0.8256, 0.7961, 0.7402, 0.5667, 0.4737, 0.3738, 0.3257, 0.2436, 0.1290, 0.0759, 0.0251, 0.0177, 0.0136, 0.0146, 0.0116, 0.0148, 0.0150, 0.0137, 0.0154, 0.0171, 0.0159, 0.0298],
[6, 0.8580, 0.7877, 0.7964, 0.7885, 0.7738, 0.7345, 0.6076, 0.4787, 0.3220, 0.2507, 0.2466, 0.2400, 0.0263, 0.0176, 0.0150, 0.0131, 0.0122, 0.0132, 0.0144, 0.0149, 0.0155, 0.0157, 0.0251, 0.0427],
[8, 0.7551, 0.7685, 0.7402, 0.7331, 0.6331, 0.6016, 0.6152, 0.6482, 0.7130, 0.6816, 0.7180, 0.8051, 0.0307, 0.0184, 0.0156, 0.0145, 0.0136, 0.0148, 0.0144, 0.0146, 0.0153, 0.0177, 0.0230, 0.0365]]
)

theta = np.arange(10,180,15)
print(theta)

mpl.rc("text", usetex=True)
nrows=3
ncols=1
fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=(3.69,5.))
#plt.rcParams.update({'font.size': 10})

plt.subplot(nrows,ncols,1)

plt.errorbar( theta, S0_v0[0,1:len(theta)+1], S0_v0[0,len(theta)+1:], label=r'$J=0$', capsize=4 )
plt.errorbar( theta, S0_v0[1,1:len(theta)+1], S0_v0[1,len(theta)+1:], label=r'$J=2$', capsize=4 )
plt.errorbar( theta, S0_v0[2,1:len(theta)+1], S0_v0[2,len(theta)+1:], label=r'$J=4$', capsize=4 )
plt.errorbar( theta, S0_v0[3,1:len(theta)+1], S0_v0[3,len(theta)+1:], label=r'$J=6$', capsize=4 )
plt.errorbar( theta, S0_v0[4,1:len(theta)+1], S0_v0[4,len(theta)+1:], label=r'$J=8$', capsize=4 )

plt.annotate(r'$\nu=0$', xy=(0.5,0.1), xycoords='axes fraction', horizontalalignment='center')

plt.legend(loc='best', numpoints=1, frameon=True, ncol=2, columnspacing=1., fontsize=8)
plt.xlim(0., 180.)
plt.xticks(np.linspace(0.,180.,5))
plt.ylim(0.0, 0.9)
#plt.ylim(1e-4, 1.)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False)
plt.ylabel('Reaction probability')
#plt.xlabel(r'$\theta$ (degrees)')

plt.subplot(nrows,ncols,2)

plt.errorbar( theta, S0_v1[0,1:len(theta)+1], S0_v1[0,len(theta)+1:], label=r'$J=0$', capsize=4 )
plt.errorbar( theta, S0_v1[1,1:len(theta)+1], S0_v1[1,len(theta)+1:], label=r'$J=2$', capsize=4 )
plt.errorbar( theta, S0_v1[2,1:len(theta)+1], S0_v1[2,len(theta)+1:], label=r'$J=4$', capsize=4 )
plt.errorbar( theta, S0_v1[3,1:len(theta)+1], S0_v1[3,len(theta)+1:], label=r'$J=6$', capsize=4 )
plt.errorbar( theta, S0_v1[4,1:len(theta)+1], S0_v1[4,len(theta)+1:], label=r'$J=8$', capsize=4 )

plt.annotate(r'$\nu=1$', xy=(0.5,0.1), xycoords='axes fraction', horizontalalignment='center')

#plt.legend(loc='best', numpoints=1, frameon=True, ncol=2, columnspacing=1.)
plt.xlim(0., 180.)
plt.xticks(np.linspace(0.,180.,5))
plt.ylim(0.0, 0.9)
#plt.ylim(1e-4, 1.)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False)
plt.ylabel('Reaction probability')
#plt.xlabel(r'$\theta$ (degrees)')

plt.subplot(nrows,ncols,3)

plt.errorbar( theta, S0_v2[0,1:len(theta)+1], S0_v2[0,len(theta)+1:], label=r'$J=0$', capsize=4 )
plt.errorbar( theta, S0_v2[1,1:len(theta)+1], S0_v2[1,len(theta)+1:], label=r'$J=2$', capsize=4 )
plt.errorbar( theta, S0_v2[2,1:len(theta)+1], S0_v2[2,len(theta)+1:], label=r'$J=4$', capsize=4 )
plt.errorbar( theta, S0_v2[3,1:len(theta)+1], S0_v2[3,len(theta)+1:], label=r'$J=6$', capsize=4 )
plt.errorbar( theta, S0_v2[4,1:len(theta)+1], S0_v2[4,len(theta)+1:], label=r'$J=8$', capsize=4 )

plt.annotate(r'$\nu=2$', xy=(0.5,0.1), xycoords='axes fraction', horizontalalignment='center')

#plt.legend(loc='best', numpoints=1, frameon=True, ncol=2, columnspacing=1.)
plt.xlim(0., 180.)
plt.xticks(np.linspace(0.,180.,5))
plt.ylim(0.0, 0.9)
#plt.ylim(1e-4, 1.)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True)
plt.ylabel('Reaction probability')
plt.xlabel(r'$\theta_\textrm{i}$ (degrees)')

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