#!/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)
#Nrows = 3
#Ncols = 1
#fig, ax = plt.subplots(nrows=Nrows, ncols=Ncols, figsize=(3.69,4.))
fig = plt.figure(figsize=(3.69, 5.))

outer = mpl.gridspec.GridSpec(2, 1, height_ratios = [3, 1])
gs1 = mpl.gridspec.GridSpecFromSubplotSpec(3, 1, subplot_spec = outer[0], hspace=0.)
gs2 = mpl.gridspec.GridSpecFromSubplotSpec(1, 1, subplot_spec = outer[1])

#Ei, ratio, standard error
NN_533_D2 = np.array(
[[0.0073, 0.4361, 0.0095, 0.0000, 0.0471],
[0.0304, 0.5457, 0.0155, 0.0000, 0.0142],
[0.1436, 0.4658, 0.0256, 0.3624, 0.0247]]
)

NN_322_D2 = np.array(
[[0.0073, 0.3134, 0.0089, 0.0000, 0.0238],
[0.0304, 0.3329, 0.0173, 0.0400, 0.0277],
[0.1436, 0.3328, 0.0193, 0.2459, 0.0175]]
)

NN_755_D2 = np.array(
[[0.0073, 0.4514, 0.0095, 0.0714, 0.0487],
[0.0304, 0.3599, 0.0175, 0.2956, 0.0362],
[0.1436, 0.4559, 0.0173, 0.4515, 0.0170]]
)

NN_433_D2 = np.array(
[[0.0073, 0.3800, 0.0030, 0.0918, 0.0104],
[0.0304, 0.3537, 0.0056, 0.3070, 0.0096],
[0.1436, 0.4034, 0.0047, 0.3740, 0.0046]]
)

NN_977_D2 = np.array(
[[0.0073, 0.4254, 0.0097, 0.3248, 0.0433],
[0.0304, 0.3715, 0.0172, 0.3697, 0.0286],
[0.1436, 0.4758, 0.0146, 0.4698, 0.0145]]
)

NN_544_D2 = np.array(
[[0.0073, 0.4187, 0.0097, 0.3019, 0.0364],
[0.0304, 0.3719, 0.0240, 0.3886, 0.0336],
[0.1436, 0.4440, 0.0131, 0.4273, 0.0130]]
)

NN_1199_D2 = np.array(
[[0.0073, 0.4582, 0.0099, 0.4619, 0.0355],
[0.0304, 0.4083, 0.0169, 0.4178, 0.0239],
[0.1436, 0.4502, 0.0128, 0.4455, 0.0128]]
)

NN_655_D2 = np.array(
[[0.0073, 0.4606, 0.0098, 0.3610, 0.0309],
[0.0304, 0.3834, 0.0170, 0.3815, 0.0236],
[0.1436, 0.4654, 0.0125, 0.4534, 0.0124]]
)

NN_131111_D2 = np.array(
[[0.0073, 0.4580, 0.0100, 0.4427, 0.0312],
[0.0304, 0.4302, 0.0168, 0.4397, 0.0219],
[0.1436, 0.4771, 0.0118, 0.4789, 0.0118]]
)

NN_766_D2 = np.array(
[[0.0073, 0.4424, 0.0101, 0.3736, 0.0297],
[0.0304, 0.4430, 0.0165, 0.4523, 0.0209],
[0.1436, 0.4560, 0.0118, 0.4472, 0.0118]]
)

NN_151313_D2 = np.array(
[[0.0073, 0.4812, 0.0103, 0.4984, 0.0283],
[0.0304, 0.4730, 0.0158, 0.4788, 0.0198],
[0.1436, 0.4870, 0.0118, 0.4871, 0.0117]]
)

NN_877_D2 = np.array(
[[0.0073, 0.4457, 0.0078, 0.4084, 0.0202],
[0.0304, 0.4482, 0.0165, 0.4806, 0.0210],
[0.1436, 0.4601, 0.0114, 0.4564, 0.0113]]
)

NN_171515_D2 = np.array(
[[0.0073, 0.4653, 0.0104, 0.4601, 0.0276],
[0.0304, 0.4664, 0.0167, 0.4822, 0.0206],
[0.1436, 0.4787, 0.0112, 0.4795, 0.0112]]
)

# 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': 1./2.,
'322': 1./3.,
'755': 2./4.,
'433': 2./5.,
'977': 3./6.,
'544': 3./7.,
'1199': 4./8.,
'655': 4./9.,
'131111': 5./10.,
'766': 5./11.,
'151313': 6./12.,
'877': 6./13.,
'171515': 7./14.
}

stat_density = []
stat_density.append( NN_step_density['533'] )
stat_density.append( NN_step_density['322'] )
stat_density.append( NN_step_density['755'] )
stat_density.append( NN_step_density['433'] )
stat_density.append( NN_step_density['977'] )
stat_density.append( NN_step_density['544'] )
stat_density.append( NN_step_density['1199'] )
stat_density.append( NN_step_density['655'] )
stat_density.append( NN_step_density['131111'] )
stat_density.append( NN_step_density['766'] )
stat_density.append( NN_step_density['151313'] )
stat_density.append( NN_step_density['877'] )
stat_density.append( NN_step_density['171515'] )
stat_ratio = []
stat_ratio.append( NN_step_ratio['533'] )
stat_ratio.append( NN_step_ratio['322'] )
stat_ratio.append( NN_step_ratio['755'] )
stat_ratio.append( NN_step_ratio['433'] )
stat_ratio.append( NN_step_ratio['977'] )
stat_ratio.append( NN_step_ratio['544'] )
stat_ratio.append( NN_step_ratio['1199'] )
stat_ratio.append( NN_step_ratio['655'] )
stat_ratio.append( NN_step_ratio['131111'] )
stat_ratio.append( NN_step_ratio['766'] )
stat_ratio.append( NN_step_ratio['151313'] )
stat_ratio.append( NN_step_ratio['877'] )
stat_ratio.append( NN_step_ratio['171515'] )

def plot_ratio( idx, title ):
	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 = []
	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='C0{:d}'.format(idx), label=title, capsize=4 )
	
	plt.plot( stat_density, stat_ratio, marker='', color='k', ls='--', label='Statistical' )
	
	plt.annotate(title, xy=(0.05, 0.2), xycoords='axes fraction', size=9)

	return

ax = plt.subplot(gs1[0])

plot_ratio(0, r'(a) $E_\textrm{i}=0.7$ kJ/mol')

#plt.legend(loc='best', numpoints=1, frameon=True, fontsize=8)
plt.xlim(0., 1.2)
plt.ylim(0.0, 0.6)
#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)')

ax = plt.subplot(gs1[1])

plot_ratio(1, r'(b) $E_\textrm{i}=2.9$ kJ/mol')

#plt.legend(loc='best', numpoints=1, frameon=True, fontsize=8)
plt.xlim(0., 1.2)
plt.ylim( 0., 0.59 )
#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(r'Ratio Terrace$_\textrm{shadow}$/Terrace$_\textrm{total}$')

ax = plt.subplot(gs1[2])

plot_ratio(2, r'(c) $E_\textrm{i}=13.9$ kJ/mol')

#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.59)
#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}$)')

ax = plt.subplot(gs2[0])

def plot_barriers( B, label, idx ):
	xaxis = np.linspace(1, len(B), len(B))
	plt.plot( xaxis[0], B[0], marker='o', color='C{:02d}'.format(idx), label=r'{:3.2f} ({:s})'.format( NN_step_density[label], label ), zorder=idx )
	plt.plot( xaxis, B, color='C{:02d}'.format(idx), zorder=idx )
	plt.scatter( xaxis[::2], B[::2], marker='o', edgecolor='C{:02d}'.format(idx), zorder=idx )
	plt.scatter( xaxis[1::2], B[1::2], marker='o', edgecolor='C{:02d}'.format(idx), facecolor='None', zorder=idx )
	

B433 = np.array([-2.175279, 1.223509, -2.20282, 3.255121, 4.017645])
B544 = np.array([-4.373374, 0.231963, -1.925594, -3.628178, -3.08126, 2.224665, 2.948515])
B655 = np.array([-5.41743, -1.288884, -3.412598, -3.138706, -3.550043, -4.083387, -3.7594, 1.672128, 2.535156])
B766 = np.array([-6.987767, -3.067319, -4.189405, -5.07005, -3.923468, -4.216894, -4.084674, -4.566505, -3.534343, 1.443053, 2.266226])
B877 = np.array([-6.323673, -2.193185, -4.658713, -4.781586, -4.542613, -4.64296, -4.458545, -4.575773, -4.328819, -4.735876, -4.199986, 1.279065, 2.16076])

idx = 0
plot_barriers( B433, '433', idx ); idx += 1
plot_barriers( B544, '544', idx ); idx += 1
plot_barriers( B655, '655', idx ); idx += 1
plot_barriers( B766, '766', idx ); idx += 1
plot_barriers( B877, '877', idx ); idx += 1

plt.annotate(r'(d)', xy=(0.05, 0.85), xycoords='axes fraction', size=9)

plt.legend(loc='best', numpoints=1, frameon=True, fontsize=7)
plt.xlim(0.5, 18.5)
plt.xticks( np.arange(1, 15, 2) )
plt.ylim(-8, 4.5)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True)
plt.ylabel(r'Barrier (kJ/mol)')
plt.xlabel(r'Terrace top site')

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