function p=newtoninterp(xx,x,a)
%Calculate values using Newton's polynomial interpolation
% x -> column vector of values along x-axis
% a -> Newton interpolating coefficients for polynomial fit to data
% xx -> column vector of interpolation points required
n=length(x);
m=length(xx);
ddd(1:m,1)=diag(eye(m));
c(1:m,1)=a(1);
for i=2:n
   ddd(1:m,i)=(xx(1:m)-x(i-1)).*ddd(1:m,i-1);
   c(1:m,i)=a(i).*ddd(1:m,i);
   p(1:m)=sum(c,2);
end
