# polynomial interpolation according to Neville's algorithm,
# i.e. "Numerical Recipes" polint() on page 109 translated into Python.
# polint() takes two lists, xa, ya, size n, and input x, and pass out
# two lists of size 1 to emulate C for y and dy. y is interpolated value
# and dy is its error. Unlike the C code, here list index starts from 0.
#
import math
def polint(xa, ya, n, x, y, dy):
# initialize c and d to be equal to ya
c = ya.copy()
d = ya.copy()
# then find an index ns which is closest to the value x in xa()
ns = 0
dif = math.fabs(x-xa[0])
for i in range(1,n):
dift = math.fabs(x-xa[i])
if (dift < dif):
ns = i
dif = dift
y[0] = ya[ns]
ns -= 1
# do double loop over column m and row i on the triangular Tableau.
# we don't try to catch the dividing by 0 error, let Python itself do it
for m in range(1,n):
for i in range(0,n-m):
ho = xa[i]-x
hp = xa[i+m]-x
w = c[i+1]-d[i]
den = ho - hp
den = w/den
d[i] = hp*den
c[i] = ho*den
# end for i
# tricking coding here in C
if(2*(ns+1) < (n-m)):
dy[0] = c[ns+1]
else:
dy[0] = d[ns]
ns -= 1
# end if else
y[0] += dy[0]
# end for m
print(c)
print(d)
# end def polint(...)
# a test run
xa = [0,1,2,4]
ya = [1,2,3,0]
x = 3.0
y = [0]
dy = [0]
n = 4
polint(xa,ya,n,x,y,dy)
print("y=",y, "dy=",dy)