# -*- coding: utf-8 -*-
"""
Created on Sat Jan 14 18:10:51 2017

@author: Mariusz
"""

# Translation to Python/Numpy of Sauer's
# Program 3.1 Newton Divided Difference Interpolation Method
# Computes coefficients of interpolating polynomial
# Input: x and y are vectors containing the x and y coordinates
#        of the n data points
# Output: coefficients c of interpolating polynomial in nested form
# Use with nest.m to evaluate interpolating polynomial

from numpy import *

def newtdd(x, y):
    n = len(x)
    v = zeros((n, n))
    for j in range(n):
        v[j,0] = y[j]			# Fill in y column of Newton triangle
    for i in range(1,n):		# For column i,
        for j in range(n-i):	# 1:n+1-i		# fill in column from top to bottom
#		print j,i," ",j+1,i-1," ", j,i-1," ",j+i," ",j
            v[j,i] = (v[j+1,i-1]-v[j,i-1])/(x[j+i]-x[j])
    c = v[0,:].copy()			# Read along top of triangle for output coefficients
    return c

'''
# test_newtdd.py
from newtdd import newtdd
from nest   import nest
from numpy  import linspace
from pylab  import plot, show
xdata = [0.0,1,2,3,3.4]
ydata = [0.0,3,12,6,2]
plot(xdata,ydata,'ro')
interpolant = newtdd(xdata,ydata)
x = linspace(-1,4,100,endpoint=True)
y = nest(interpolant,x,xdata)
plot(x,y,'b')
show()
'''