% Programming and Numerical Methods
% MSc in Computational and Applied Physics

% Code 7: An non-equispaced interpolation
% Zeros of Chebyshev polynomials
% A. Meseguer (2007)

clear all ; close all 
% As in Code 1, we first set up the Lagrange poly.

 N = 32 ;                   % Number of nodes in = 0,1,...,N
 a = -5 ; b = 5 ;           % Domain: x in [a,b]
 xn = (a+b)/2 - 0.5*(b-a)*cos(pi*[0:N]'/N); % Spatial grid: [x_0, x_1, ... , x_N]' - N x 1  

 NN = 200 ; xx = linspace(a,b,NN)' ;  % A finer grid for evaluation - 100 x 1

 CP = ones(NN,N+1) ;      % CP is full of 1's 
 for ipol = 0:N
   fnorm = 1 ;
   for jprod = 0:N
    if jprod ~= ipol ; fnorm = fnorm*(xn(ipol+1)-xn(jprod+1)) ; end
   end
   
   for jprod = 0:N
    if jprod ~= ipol ; CP(:,ipol+1) = CP(:,ipol+1).*(xx - xn(jprod+1)) ; end
   end
   CP(:,ipol+1) = CP(:,ipol+1)/fnorm ;
 end

% Now we define the function to be interpolated
  f = exp(-xn.^2); % Original function at the interpolation points
  ff = exp(-xx.^2) ; % The same, but on a finer grid for plot quality

% First: we plot the function in a finer grid
  plot(xx,ff,'-b') ; axis([-5 5 -1.5 1.5]) ; grid on ; hold on

% Second: we represent the points where the function is going to be
% interpolated
  plot(xn,f,'or','markerfacecolor',[1 0 0]) ; grid on ;

% Third: we plot interpolating polynomial
  PIN = CP*f ;
  plot(xx,PIN,'-m')

% Lets have a look at the error
  figure(2) 
  semilogy(xx,abs(PIN-ff),'ok','markerfacecolor',[0 0 0])
  title('ERROR')

