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

# E, top, hcp, fcc, bridge, total
n1100K = np.array(
[[74.0, 6, 2, 15, 23, 500.],
[87.8, 6, 9, 31, 46, 500.],
[102.9, 21, 14, 38, 73, 500.],
[119.5, 11, 24, 51, 86, 500.]]
)

n475K = np.array(
[[102.9, 28, 36, 68, 132, 1000.],
[119.5, 29, 41, 124, 194, 1000.]]
)

stick_1100K = np.array(
[[74.0, 0.046, 0.009, 0.196, 0.018],
[87.8, 0.092, 0.013, 0.192, 0.018],
[102.9, 0.146, 0.016, 0.220, 0.019],
[119.5, 0.172, 0.017, 0.224, 0.019]]
)

stick_475K = np.array(
[[102.9, 0.132, 0.011, 0.209, 0.013],
[119.5, 0.194, 0.013, 0.244, 0.014]]
)

def confidence(P,T):
        'compute the 1 sigma confidence given a probability P and the total number of events T'
	if P == 0:
		sigma = 1 - 0.68**(1./T)
	else:
		sigma = np.sqrt( P * (1-P) / T )

        return sigma

def site( nsites ):
	psites = np.zeros((len(nsites), 7))

	for i in range(len(nsites)):
		psites[i][0] = nsites[i][0]

		psites[i][1] = nsites[i][1] / nsites[i][4]
		psites[i][2] = confidence( psites[i][1], nsites[i][4] )

		psites[i][3] = nsites[i][2] / nsites[i][4]
		psites[i][4] = confidence( psites[i][3], nsites[i][4] )

		psites[i][5] = nsites[i][3] / nsites[i][4]
		psites[i][6] = confidence( psites[i][5], nsites[i][4] )

	return psites

def prob( nsites ):
	psites = np.zeros((len(nsites), 7))

	for i in range(len(nsites)):
		psites[i][0] = nsites[i][0]

		psites[i][1] = nsites[i][1] / (nsites[i][5])
		psites[i][2] = confidence( psites[i][1], nsites[i][5] )

		psites[i][3] = nsites[i][2] / (nsites[i][5])
		psites[i][4] = confidence( psites[i][3], nsites[i][5])

		psites[i][5] = nsites[i][3] / (nsites[i][5])
		psites[i][6] = confidence( psites[i][5], nsites[i][5])

	return psites

site1100K = site( n1100K )
site475K = site( n475K )
p1100K = prob( n1100K )
p475K = prob( n475K )

mpl.rc("text", usetex=True)
fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(3.69,4.))

plt.subplot(2,1,1)

plt.errorbar( site1100K[:,0], site1100K[:,1], yerr=site1100K[:,2], marker='o', label='Top', color='b')
plt.errorbar( site1100K[:,0], site1100K[:,3], yerr=site1100K[:,4], marker='d', label='Hollow', color='r')
plt.errorbar( site1100K[:,0], site1100K[:,5], yerr=site1100K[:,6], marker='v', label='Bridge', color='g')

plt.errorbar( site475K[:,0], site475K[:,1], yerr=site475K[:,2], marker='o', color='b', fillstyle='none', linestyle='--')
plt.errorbar( site475K[:,0], site475K[:,3], yerr=site475K[:,4], marker='d', color='r', fillstyle='none', linestyle='--')
plt.errorbar( site475K[:,0], site475K[:,5], yerr=site475K[:,6], marker='v', color='g', fillstyle='none', linestyle='--')

plt.plot([0., 400.], [0.25,0.25], linestyle=(0, (1, 4)), color='b', label='')
plt.plot([0., 400.], [0.25,0.25], linestyle=(0, (1, 9)), color='r', label='')
plt.plot([0., 400.], [0.5,0.5], linestyle=':', color='g', label='')

#plt.legend(loc='upper left', numpoints=1, handletextpad=0.5, borderaxespad=0.2, frameon=True)
plt.xlim(30.,125.)
plt.ylim(0.,0.8)
plt.tick_params(labelbottom='off', length=6, width=1)
#plt.xlabel('Incidence energy (kJ/mol)')
plt.ylabel('Probability of closest site')
plt.annotate('(a)', xy=(40, 0.65))

ax2 = plt.subplot(2,1,2)

#plt.errorbar( stick_1100K[:,0], stick_1100K[:,1], yerr=stick_1100K[:,2], marker='s', label='Total', color='k')
plt.errorbar( p1100K[:,0], p1100K[:,1], yerr=p1100K[:,2], marker='o', label='Top', color='b')
plt.errorbar( p1100K[:,0], p1100K[:,3], yerr=p1100K[:,4], marker='d', label='Hollow', color='r')
plt.errorbar( p1100K[:,0], p1100K[:,5], yerr=p1100K[:,6], marker='v', label='Bridge', color='g')

#plt.errorbar( stick_475K[:,0], stick_475K[:,1], yerr=stick_475K[:,2], marker='s', color='k', fillstyle='none')
plt.errorbar( p475K[:,0], p475K[:,1], yerr=p475K[:,2], marker='o', color='b', fillstyle='none', linestyle='--')
plt.errorbar( p475K[:,0], p475K[:,3], yerr=p475K[:,4], marker='d', color='r', fillstyle='none', linestyle='--')
plt.errorbar( p475K[:,0], p475K[:,5], yerr=p475K[:,6], marker='v', color='g', fillstyle='none', linestyle='--')

plt.legend(loc='upper left', numpoints=1, handletextpad=0.5, borderaxespad=0.2, frameon=False)
plt.xlim(30.,125.)
plt.ylim(0.,0.15)
plt.tick_params(labelleft='off', labelright='on', length=6, width=1)
ax2.yaxis.set_label_position("right")
plt.xlabel('Incidence energy (kJ/mol)')
plt.ylabel('Site specific sticking probability')
plt.annotate('(b)', xy=(100, 0.13))

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