% --------------------------------------------------------
% RBF Interpolation with gaussians, monomials and TPS in 1D
%
% Test function: f(x)=sin(x^2/2) x in [0,pi]
%
% Author: Stefano De Marchi
%             March 19, 2013
%----------------------------------------------------------

clear all;
close all;

x=linspace(0,pi,10);
%f=[.2 3.4 -1 0 .6];
f=sin(x.^2/2);
d=input('give me the polynomial degree : '); 
m=length(x);

if d>m-1
    error('degree too big !');
end

% Define the Gaussian RBF, monomials, TPS and shape parameter
  rbf  = inline('exp(-(e*r).^2)');
  rbf1 = inline('r.^d');
  rbf2 = inline('r.^2 .* log(r+e)');
  
  
  %ep=input('give me the shape parameter = ');
  ep = 2;
 
  neval = 200; 
  grid = linspace(0,pi,neval);
 
  % Evaluate the test function at the data points
  rhs = f';
  % Compute distance matrix between the data sites and centers
  DM_data = DistanceMatrix(x',x');
  % Compute interpolation matrix
  IM = feval(rbf,ep,DM_data);
  IM1 = feval(rbf1,d,DM_data);
  IM2 = feval(rbf2,ep,DM_data);
  % Compute distance matrix between evaluation points and centers
  DM_eval = DistanceMatrix(grid',x');
  % Compute evaluation matrix
  EM = feval(rbf,ep,DM_eval);
  EM1 = feval(rbf1,d,DM_eval);
  EM2 = feval(rbf2,ep,DM_eval);
  % Compute RBF interpolant
  % (evaluation matrix * solution of interpolation system)
  Pf = EM * (IM\rhs);
  Pf1 = EM1 * (IM1\rhs);
  Pf2 = EM2 * (IM2\rhs);

xx=grid;
ff=sin(xx.^2/2);

if (d==3)
    sp=spline(x,f,xx);
    plot(xx,ff,'-k',xx,Pf,'-b',xx,Pf1,'-r',xx,Pf2,'g-',xx,sp,'c-');
    legend('funzione','G', 'MN','TPS','spline');
    stringa=['polynomial degree ', num2str(d)];
    title(stringa);
else
    plot(xx,ff,'-k',xx,Pf,'-b',xx,Pf1,'r-',xx,Pf2,'g-');
    legend('funzione','G','interp. MN','TPS');
    stringa=['polynomial degree ', num2str(d)];
    title(stringa);
end