# Integration of tabular data using trapezoidal 
# and Simpson's rule

import numpy as np
from scipy import integrate

# 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.7614402697

# Integrate table values using trapz function
# Note the order of the inputs is backwards from
# the MATLAB trapz function
t = np.trapz(y,x)
print('t = ',t)
ret = (t-exact)/exact
print('ret = ',ret)

# If table has odd number of values you can also
# do
# s = integrate.simpson(y,x)
# print('s = ',s)

# If the table has an even number of values, you
# can still use Simpson's rule by specifying one
# of three ways to handle the situation

# first option uses Simpson's for first points
# and a trapezoidal rule on last panel
s1 = integrate.simpson(y,x,even='first')
print('s1 = ',s1)
res1 = (s1-exact)/exact
print('ret = ',res1)

# last option uses Simpson's for last points
# and a trapezoidal rule on first panel
s2 = integrate.simpson(y,x,even='last')
print('s2 = ',s2)
res2 = (s2-exact)/exact
print('ret = ',res2)

# ave option does both of the above and averages the two
# results.
s3 = integrate.simpson(y,x,even='avg')
print('s3 = ',s3)
res3 = (s3-exact)/exact
print('ret = ',res3)