#!/usr/bin/env python

######################################################### 
import 	matplotlib.pyplot as plt			#	
import 	matplotlib        as mpl			#
import 	numpy             as np				#
from   	math              import *			#

# set up plot parameters
mpl.rc("font", size=10)
fig = plt.figure(figsize=(3.69,3))

def plotfunction( plot, filename ):
	MEP_Z = []
	MEP_r = []
	MEP_E = []

	lines = open(filename).readlines()
	for i in range(1, len(lines)):
		MEP_r.append(float(  lines[i].split()[0]  ))
		MEP_Z.append(float(  lines[i].split()[1]  ))
		MEP_E.append(float(  lines[i].split()[2]  )) 
	pos_max = np.argmax(MEP_E)

	s = [ 0. ]
	for i in range( 1, len( MEP_E ) ):
		dZ = MEP_Z[i] - MEP_Z[i-1]
		dr = MEP_r[i] - MEP_r[i-1]
		if i == 1:
			s.append( sqrt( dZ**2 + dr**2 ) )
		else:
			s.append( s[-1] + sqrt( dZ**2 + dr**2 ) )
	ds = s[pos_max]
	for i in range( len(s) ):
		s[i] -= ds

	if plot == 1:
		plt.scatter(s, MEP_E, color = 'blue', zorder=1, label='Cu(111)')
		plt.scatter(s[pos_max], MEP_E[pos_max], color = 'black', marker='s', s=35, zorder=2)
	elif plot == 2:
		plt.scatter(s, MEP_E, color = 'red', zorder=1, label='Ni(111)')
		plt.scatter(s[pos_max], MEP_E[pos_max], color = 'black', marker='s', s=35, zorder=2)
		plt.xlim(-1.5, 0.6)
		plt.ylim(0., 200.)
		plt.xlabel( r's ($\AA$)' )
		plt.ylabel( r'E (kJ/mol)' )
#		plt.annotate('(d)', xy=(1.3, 2.8), size=16)
		plt.legend( scatterpoints=1, loc='upper left' )

	return

plotfunction( 1, 'energy_nn_mep.dat' )
plotfunction( 2, 'energy_ni_mep.dat' )

plt.tight_layout()

plt.savefig('mep_ni.pdf')
#plt.show()