# Integration of tabular data using cubic spline
# integration

import numpy as np
from scipy import interpolate 
import matplotlib.pyplot as plt

# Use the loadtxt function to load the data into
# an np array.  This needs to be transposed so that
# the extraction of columns can be done.
A = np.loadtxt('hw31.dat')
A = np.transpose(A)

# Extract first column into x and second column into
# y
x = A[0]
y = A[1]
exact = 0.76144026974 

# Compute the cubic spline structure
#  k = 1 = linear spline
#  k = 2 = quadratic spline
#  k = 3 = cubic spline (not-a-knot)
S = interpolate.splrep(x,y,k=3)

# Grow to a new table with more points; plot both
# of these
n = len(x)
xo = np.linspace(x[0],x[n-1],51)
yo = interpolate.splev(xo,S,der=0)
plt.figure(1)
plt.plot(x,y,'b*')
plt.plot(xo,yo,'ro')

# Integrate the spline.  s1 is a tuple here also
s1 = interpolate.splint(x[0],x[n-1],S)
re = (s1-exact)/exact
print('Relative error in integral',re)