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

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

#Ei, ratio, standard error
NN_533_D2 = np.array(
[[0.0073, 0.6046, 0.0059, 0.9988, 0.0004],
[0.0304, 0.7371, 0.0070, 0.9931, 0.0013],
[0.1436, 0.8625, 0.0065, 0.8632, 0.0065]]
)

NN_322_D2 = np.array(
[[0.0073, 0.5424, 0.0065, 0.9973, 0.0007],
[0.0304, 0.7470, 0.0080, 0.9830, 0.0024],
[0.1436, 0.7659, 0.0084, 0.7616, 0.0085]]
)

NN_755_D2 = np.array(
[[0.0073, 0.4791, 0.0069, 0.9947, 0.0010],
[0.0304, 0.6914, 0.0094, 0.9348, 0.0050],
[0.1436, 0.6570, 0.0097, 0.6454, 0.0097]]
)

NN_433_D2 = np.array(
[[0.0073, 0.4496, 0.0023, 0.9839, 0.0006],
[0.0304, 0.6671, 0.0032, 0.8950, 0.0021],
[0.1436, 0.5507, 0.0032, 0.5455, 0.0032]]
)

NN_977_D2 = np.array(
[[0.0073, 0.4276, 0.0074, 0.9741, 0.0024],
[0.0304, 0.6183, 0.0107, 0.8621, 0.0076],
[0.1436, 0.5134, 0.0102, 0.5072, 0.0102]]
)

NN_544_D2 = np.array(
[[0.0073, 0.3744, 0.0076, 0.9612, 0.0030],
[0.0304, 0.5592, 0.0164, 0.7709, 0.0138],
[0.1436, 0.4299, 0.0099, 0.4236, 0.0099]]
)

NN_1199_D2 = np.array(
[[0.0073, 0.3501, 0.0076, 0.9497, 0.0035],
[0.0304, 0.5313, 0.0118, 0.7637, 0.0100],
[0.1436, 0.3685, 0.0098, 0.3685, 0.0098]]
)

NN_655_D2 = np.array(
[[0.0073, 0.3287, 0.0076, 0.9371, 0.0039],
[0.0304, 0.5205, 0.0121, 0.7529, 0.0104],
[0.1436, 0.3433, 0.0096, 0.3413, 0.0096]]
)

NN_131111_D2 = np.array(
[[0.0073, 0.3113, 0.0077, 0.9299, 0.0042],
[0.0304, 0.4704, 0.0123, 0.6860, 0.0115],
[0.1436, 0.2789, 0.0090, 0.2757, 0.0090]]
)

NN_766_D2 = np.array(
[[0.0073, 0.2988, 0.0078, 0.9224, 0.0046],
[0.0304, 0.4370, 0.0124, 0.6471, 0.0119],
[0.1436, 0.2834, 0.0090, 0.2809, 0.0090]]
)

NN_151313_D2 = np.array(
[[0.0073, 0.2874, 0.0079, 0.9055, 0.0051],
[0.0304, 0.3947, 0.0120, 0.6144, 0.0120],
[0.1436, 0.2669, 0.0089, 0.2637, 0.0089]]
)

NN_877_D2 = np.array(
[[0.0073, 0.2917, 0.0060, 0.8958, 0.0040],
[0.0304, 0.3975, 0.0126, 0.6231, 0.0125],
[0.1436, 0.2329, 0.0085, 0.2285, 0.0084]]
)

NN_171515_D2 = np.array(
[[0.0073, 0.2633, 0.0079, 0.8958, 0.0055],
[0.0304, 0.3881, 0.0127, 0.5969, 0.0128],
[0.1436, 0.2176, 0.0082, 0.2157, 0.0081]]
)

# In terms of unit cell vector parallel to the terrace and in nm^-1
# Please check the non-rectangular unit cells, there might be an error (or not)
NN_step_density = {
'211': 1.33925980477299,
'533': 1.10245499657604,
'322': 0.867222846526904,
'755': 0.717738278424414,
'433': 0.609594900345529,
'977': 0.531226638707765,
'544': 0.469525140800388,
'1199': 0.42152023770247,
'655': 0.381732224349042,
'131111': 0.349345915692593,
'766': 0.321582789240287,
'151313': 0.298239412749505,
'877': 0.277790623773609,
'171515': 0.260184400328028
}

NN_step_ratio = {
# '211': ,
'533': 0.50400,
'322': 0.39950,
'755': 0.33065,
'433': 0.28160,
'977': 0.24495,
'544': 0.21660,
'1199': 0.19400,
'655': 0.17565,
'131111': 0.16050,
'766': 0.14770,
'151313': 0.13670,
'877': 0.12740,
'171515': 0.11910
}

def plot_ratio( idx ):
	density = []
	ratio = []
	density.append( NN_step_density['533'] )
	density.append( NN_step_density['322'] )
	density.append( NN_step_density['755'] )
	density.append( NN_step_density['433'] )
	density.append( NN_step_density['977'] )
	density.append( NN_step_density['544'] )
	density.append( NN_step_density['1199'] )
	density.append( NN_step_density['655'] )
	density.append( NN_step_density['131111'] )
	density.append( NN_step_density['766'] )
	density.append( NN_step_density['151313'] )
	density.append( NN_step_density['877'] )
	density.append( NN_step_density['171515'] )
	ratio.append( [ NN_533_D2[idx,1], NN_533_D2[idx,2] ] )
	ratio.append( [ NN_322_D2[idx,1], NN_322_D2[idx,2] ] )
	ratio.append( [ NN_755_D2[idx,1], NN_755_D2[idx,2] ] )
	ratio.append( [ NN_433_D2[idx,1], NN_433_D2[idx,2] ] )
	ratio.append( [ NN_977_D2[idx,1], NN_977_D2[idx,2] ] )
	ratio.append( [ NN_544_D2[idx,1], NN_544_D2[idx,2] ] )
	ratio.append( [ NN_1199_D2[idx,1], NN_1199_D2[idx,2] ] )
	ratio.append( [ NN_655_D2[idx,1], NN_655_D2[idx,2] ] )
	ratio.append( [ NN_131111_D2[idx,1], NN_131111_D2[idx,2] ] )
	ratio.append( [ NN_766_D2[idx,1], NN_766_D2[idx,2] ] )
	ratio.append( [ NN_151313_D2[idx,1], NN_151313_D2[idx,2] ] )
	ratio.append( [ NN_877_D2[idx,1], NN_877_D2[idx,2] ] )
	ratio.append( [ NN_171515_D2[idx,1], NN_171515_D2[idx,2] ] )
	ratio = np.array( ratio )

	plt.errorbar( density, ratio[:,0], yerr=ratio[:,1], marker='o', color='C00', label='Initial', capsize=4 )

	ratio = []
	ratio.append( [ NN_533_D2[idx,3], NN_533_D2[idx,4] ] )
	ratio.append( [ NN_322_D2[idx,3], NN_322_D2[idx,4] ] )
	ratio.append( [ NN_755_D2[idx,3], NN_755_D2[idx,4] ] )
	ratio.append( [ NN_433_D2[idx,3], NN_433_D2[idx,4] ] )
	ratio.append( [ NN_977_D2[idx,3], NN_977_D2[idx,4] ] )
	ratio.append( [ NN_544_D2[idx,3], NN_544_D2[idx,4] ] )
	ratio.append( [ NN_1199_D2[idx,3], NN_1199_D2[idx,4] ] )
	ratio.append( [ NN_655_D2[idx,3], NN_655_D2[idx,4] ] )
	ratio.append( [ NN_131111_D2[idx,3], NN_131111_D2[idx,4] ] )
	ratio.append( [ NN_766_D2[idx,3], NN_766_D2[idx,4] ] )
	ratio.append( [ NN_151313_D2[idx,3], NN_151313_D2[idx,4] ] )
	ratio.append( [ NN_877_D2[idx,3], NN_877_D2[idx,4] ] )
	ratio.append( [ NN_171515_D2[idx,3], NN_171515_D2[idx,4] ] )
	ratio = np.array( ratio )
	
	plt.errorbar( density, ratio[:,0], yerr=ratio[:,1], marker='o', color='C01', label='Reaction', capsize=4 )
	
	ratio = []
	ratio.append( NN_step_ratio['533'] )
	ratio.append( NN_step_ratio['322'] )
	ratio.append( NN_step_ratio['755'] )
	ratio.append( NN_step_ratio['433'] )
	ratio.append( NN_step_ratio['977'] )
	ratio.append( NN_step_ratio['544'] )
	ratio.append( NN_step_ratio['1199'] )
	ratio.append( NN_step_ratio['655'] )
	ratio.append( NN_step_ratio['131111'] )
	ratio.append( NN_step_ratio['766'] )
	ratio.append( NN_step_ratio['151313'] )
	ratio.append( NN_step_ratio['877'] )
	ratio.append( NN_step_ratio['171515'] )
	
	plt.plot( density, ratio, marker='', color='k', ls='--', label='Statistical' )
	#plt.plot( density, np.array(ratio)*2, marker='', color='k', ls=':', label='Statistical' )

	return

plt.subplot(Nrows, Ncols, 1)

plot_ratio(0)

plt.annotate(r'$E_\textrm{i}=0.7$ kJ/mol', xy=(0.6, 0.1), xycoords='axes fraction', size=8)

#plt.legend(loc='best', numpoints=1, frameon=True, fontsize=8)
plt.xlim(0., 1.2)
plt.ylim(0.0, 1.)
#plt.ylim(1e-7, 1.)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False)
#plt.ylabel('Ratio step/(step + terrace)')

plt.subplot(Nrows, Ncols, 2)

plot_ratio(1)

plt.annotate(r'$E_\textrm{i}=2.9$ kJ/mol', xy=(0.6, 0.1), xycoords='axes fraction', size=8)

#plt.legend(loc='best', numpoints=1, frameon=True, fontsize=8)
plt.xlim(0., 1.2)
plt.ylim( 0., 0.99 )
#plt.ylim(0.0, 0.8)
#plt.ylim(1e-7, 1.)
#plt.yscale('log, labelbottom=False')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False)
plt.ylabel('Ratio step/(step + terrace)')

plt.subplot(Nrows, Ncols, 3)

plot_ratio(2)

plt.annotate(r'$E_\textrm{i}=13.9$ kJ/mol', xy=(0.6, 0.1), xycoords='axes fraction', size=8)

plt.legend(loc='best', numpoints=1, frameon=True, fontsize=8)
plt.xlim(0., 1.2)
plt.ylim( 0., plt.ylim()[1] )
plt.ylim(0.0, 0.99)
#plt.ylim(1e-7, 1.)
#plt.yscale('log')
plt.tick_params(length=6, width=1, direction='in', top=True, right=True)
#plt.ylabel('Ratio step/(step + terrace)')
plt.xlabel(r'Step density (nm$^{-1}$)')

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