#!/usr/bin/env python
import sys, commands
from os import path, getcwd
from os import chdir, system, popen
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from scipy.interpolate import spline

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(3.37,3))
mpl.rc("text", usetex=True)

folders = [ '0.94', '1.18', '1.29', '1.55', '1.80', '2.12', '2.56' ]

dxy_reac = []
dxy_scat = []
dxy_trap = []

here=getcwd()
def collectxy(folder):
	chdir(folder)

	last  = int( popen('ls | egrep "^[0-9].....$" | wc -l').read() )

	here2=getcwd()				#workdir
	for i in range(1,(last+1)):
		traj= '%.6d' % i		# define working directory
		newdir=path.join(here2, traj)	# updating mypath

		if not path.exists(newdir): 
			print "folder", newdir, "not found. exiting..."
			quit()

		chdir(newdir)

		if path.exists('PostAnalysis.dat') and path.exists('outcome'):
			outcome = open('outcome', 'r').readline()
			PostDat = open('PostAnalysis.dat', 'r').readlines()
			if "SCATTERING" in outcome:
				dxy_scat.append( float( PostDat[4].split()[14] ) )
			elif "TRAPPING" in outcome:
				dxy_trap.append( float( PostDat[4].split()[14] ) )
			elif "REACTION" in outcome:
				dxy_reac.append( float( PostDat[11].split()[0] ) )

	chdir(here)

for i in range(len(folders)):
	collectxy(folders[i])

n, bins = np.histogram( dxy_reac, bins=25, density=True )
x = []
x.append( bins[0] - (bins[1]-bins[0])/2 )
for i in range(len(bins)-1):
	x.append( (bins[i]+bins[i+1])/2. )
x.append( bins[len(bins)-1] + (bins[len(bins)-1]-bins[len(bins)-2])/2 )
n = np.insert(n, 0, 0)
n = np.append(n,0)

xnew = np.linspace(0.,0.5,100.)
n_smooth = spline(x,n,xnew)
plt.plot(xnew,n_smooth, label='Reacted', c='b')

n, bins = np.histogram( dxy_scat, bins=25, density=True )
x = []
x.append( bins[0] - (bins[1]-bins[0])/2 )
for i in range(len(bins)-1):
	x.append( (bins[i]+bins[i+1])/2. )
x.append( bins[len(bins)-1] + (bins[len(bins)-1]-bins[len(bins)-2])/2 )
n = np.insert(n, 0, 0)
n = np.append(n,0)

xnew = np.linspace(0.,0.5,100.)
n_smooth = spline(x,n,xnew)
plt.plot(xnew,n_smooth, label='Scattered', c='r')

n, bins = np.histogram( dxy_trap, bins=15, density=True )
x = []
x.append( bins[0] - (bins[1]-bins[0])/2 )
for i in range(len(bins)-1):
	x.append( (bins[i]+bins[i+1])/2. )
x.append( bins[len(bins)-1] + (bins[len(bins)-1]-bins[len(bins)-2])/2 )
n = np.insert(n, 0, 0)
n = np.append(n,0)

xnew = np.linspace(0.,0.5,100.)
n_smooth = spline(x,n,xnew)
#plt.plot(xnew,n_smooth, label='Trapped', c='g')

y_limits = plt.ylim()
plt.ylim(0., y_limits[1])

plt.tick_params(labelleft='off')
plt.legend(loc='best', frameon=False)
plt.ylabel('Distribution (arb.)')
plt.xlabel(r'Steering in xy (\r{A})')

plt.tight_layout()
plt.subplots_adjust(wspace=0, hspace=0)
plt.savefig('xysteering.pdf')


print np.mean(dxy_scat)
print np.std(dxy_scat) / (len(dxy_scat)-1)**0.5
print np.mean(dxy_trap)
print np.std(dxy_trap) / (len(dxy_trap)-1)**0.5
print np.mean(dxy_reac)
print np.std(dxy_reac) / (len(dxy_reac)-1)**0.5

plt.show()