#!/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()