clear all

for n=3:15
close all

% Chebyshev points
xc=cos([0:n-1]*pi/(n-1));
[xC,yC]=ndgrid(xc);
xyC=[xC(:),yC(:)];

% equispaced points
xe=linspace(-1,1,n);
[xE,yE]=ndgrid(xe);
xyE=[xE(:),yE(:)];

% Halton points
H=2*haltonseq(n^2,2)-1;

% define the RBF
rbf = inline('exp(-(e*r).^2)');
%rbf = inline('r.^4.*(5*spones(r)-4*r)','e','r'); 
%ep=0.7*64;
ep=3;

dsites=xyE;
neval=n^2;
M=neval^2;
ctrs = dsites;   % centers coincide with data sites
grid = linspace(-1,1,neval); 
[xe,ye] = meshgrid(grid);
epoints = [xe(:) ye(:)];
  
% Compute distance matrix between the data sites and centers
  DM_data = DistanceMatrix(dsites,ctrs);
% Compute distance matrix between evaluation points and centers
  DM_eval = DistanceMatrix(epoints,ctrs);
  IM = feval(rbf,ep,DM_data);
  EM = feval(rbf,ep,DM_eval);

% Compute power function at evaluation points
     invIM = pinv(IM); 
     phi0 = rbf(ep,0);
     for j=1:M
        powfun(j) = real(sqrt(phi0-(invIM*EM(j,:)')'*EM(j,:)'));
     end
     
     PF=reshape(powfun,neval,neval);
     
     hh=surf(xe,ye,PF);
     set(hh,'FaceColor','interp','EdgeColor','none')
     % Compute max. norm of power function on evaluation grid
     maxPF = max(powfun);

     pause(2)
end
 % fprintf('Smallest maximum norm: %e\n', min(maxPF))
 % fprintf('at epsilon = %f\n',ep(find(maxPF==min(maxPF))))
 % Plot power function norm
 % figure; 
 % semilogy(ep,maxPF,'b');
