% interp_demo.m - demonstrate polynomial interpolation graphically

function interp_demo(F)
if ( nargin == 0 )
    F = @(x) 1-.9*abs(x);
end
x = [];
xx = linspace(-1,1,201)';
LW = 'linewidth'; lw = 3;
MS = 'markersize'; ms = 32;
for n = -1:100
  hold off, shg
  plot(xx,F(xx),'-k',LW,lw)               % plot function F
  hold on, axis([-1 1 -1 2])
  if n>-1
    plot(x,F(x),'.b',MS,ms)               % plot interpolation points
    plot(x,0*x,'.k',MS,ms)                 % plot x values alone
    yy = bary(x,F(x),xx);                 % compute barycentric interpolant
    plot(xx,yy,'-b',LW,lw)                % plot it
  end
  [gx,gy,button] = ginput(1);             % input new interpolation point
  if button==3, break, end                % if right button, stop
  x = unique([x; gx]);  
end

function yy = bary(x,y,xx)  % compute barycentric interpolant at points xx

% Compute barycentric weights:
  np = length(x);
  w = zeros(np,1);
  for j = 1:np
    k = 1:np; k(j) = [];
    w(j) = 1./prod(x(j)-x(k));
  end

% Apply barycentric formula:
  numer = zeros(size(xx));
  denom = zeros(size(xx));
  for i = 1:np
    xdiff = xx-x(i);
    numer = numer + w(i)*y(i)./xdiff;
    denom = denom + w(i)./xdiff;
  end
  yy = numer./denom;
