% RBFInterpolation2DCSRBF
% Script that performs 2D RBF interpolation with CSRBFs
% Calls on: DistanceMatrix
% Define a Wendland CSRBF and shape parameter
 
% Wendland C^2
rbf = inline('r.^4.*(5*spones(r)-4*r)','e','r'); 

% Wendland C^4
% rbf = inline('r.^6.*(56*spones(r)-88*r+35*r.^2)','e','r'); 

  ep=0.7; % non-stationary
  % Define Franke's function as testfunction
  f1 = '0.75*exp(-((9*x-2).^2+(9*y-2).^2)/4)';
  f2 = '0.75*exp(-((9*x+1).^2/49+(9*y+1).^2/10))';
  f3 = '0.5*exp(-((9*x-7).^2+(9*y-3).^2)/4)';
  f4 = '0.2*exp(-((9*x-4).^2+(9*y-7).^2))';
  testfunction = inline(strcat(f1,'+',f2,'+',f3,'-',f4));
  %N = 1089; gridtype = 'h';
  
  k=[3, 5, 9, 17,  33, 65];
  
  for n=1:length(k)
  close all
  % equispaced points
  xe=linspace(0,1,k(n));
  [xE,yE]=ndgrid(xe);
  xyE=[xE(:),yE(:)];
  
  %N=n^2;
  %gridtype='u';
  %Load data points
  %name = sprintf('Data2D_%d%s',N,gridtype); load(name)
  
  dsites=xyE;
  ctrs = dsites;
  neval = 40; grid = linspace(0,1,neval);
  [xe,ye] = meshgrid(grid); epoints = [xe(:) ye(:)];
  % Evaluate the test function at the data points
  rhs = testfunction(dsites(:,1),dsites(:,2));
  % Compute distance matrix between the data sites and centers
  DM_data = DistanceMatrixCSRBF_new(dsites,ctrs,ep);
  % Compute interpolation matrix
  IM = feval(rbf,ep,DM_data);
  % Compute distance matrix between evaluation points and centers
  DM_eval = DistanceMatrixCSRBF_new(epoints,ctrs,ep);
  % Compute evaluation matrix
  EM = feval(rbf,ep,DM_eval);
  
  % Compute RBF interpolant
  % (evaluation matrix * solution of interpolation system)
  sol=pcg(IM, rhs, 1.e-8, 200);
  %Pf = EM * (IM\rhs);
  Pf=EM*sol;
  
  
  % Compute exact solution, i.e.,
  % evaluate test function on evaluation points
  exact = testfunction(epoints(:,1),epoints(:,2));
  % Compute errors on evaluation grid
  maxerr = norm(Pf-exact,inf);
  rms_err = norm(Pf-exact)/neval;
  fprintf('RMS error:     %e\n', rms_err)
  fprintf('Maximum error: %e\n', maxerr)
  % Plot interpolant
  fview = [160,20]; % for Franke's function
  caption = strcat('RBF interpolant',...
      ' false colored by maximum error.');
  PlotSurf(xe,ye,Pf,neval,exact,maxerr,fview,caption);
  % Plot maximum error
  caption = 'Maximum error for RBF interpolant.';
  PlotError2D(xe,ye,Pf,exact,maxerr,neval,fview,caption);
  pause(2)
  
  %ep=0.7*2; %stationary 
  
  end