#!/usr/bin/python

import sys
from   numpy import linspace,exp
import matplotlib.pyplot as plt
from   scipy.interpolate import UnivariateSpline
from   scipy.interpolate import interp1d
from   scipy.integrate   import simps


def splineInt( x, y ):
    
   ##fit a spline without smoothing
   s  =  UnivariateSpline(x, y, s=0, k=3)
   
   ##get the integral
   dA = s.integral(0,1) 

   #dA = simps(y, x)

   return dA

f1 = open( 'C7eq/gF_all.dat', 'r')
f2 = open( 'alphaR/gF_all.dat', 'r')
f3 = open( 'C7ax/gF_all.dat', 'r')
f4 = open( 'alphaL/gF_all.dat', 'r')

slist = []
for dset in [[f1, 'r', 'C7eq'], [f2,'g','alphaR'], [f3,'b','C7ax'], [f4,'k', 'alphaL']]:

    x = []
    y = []
    
    f=dset[0]
    color=dset[1]
    lab=dset[2]
    
    for line in f:
        asList = line.split()

        ##single integration in  two-column format
        if  len(asList) == 2:
            x.append(float(line.split()[0]))
            y.append(float(line.split()[1]))

    
    dA = splineInt( x, y )
    s  =  UnivariateSpline(x, y, s=0, k=3)
    slist.append(s)
    
    print "dA %s based on spline fit:   %.6f" % (lab, dA) 
    if lab == "C7eq":
        ref = dA
    print "ddA:   %.6f\n\n" % (ref-dA) 


    xs = linspace(0,1,1000)
    ys = s(xs)

    plt.plot(x, y, color+'.')
    plt.plot(xs, ys, color+'--', label=lab)

plt.xlabel('$\lambda$', fontsize=18)
plt.ylabel('Generalised Force', fontsize=18)
plt.legend()
plt.savefig('genForces_splines.pdf')
plt.show()