clear all;
close all;
% Program for plotting the power function
% for piecewise interpolation on a triangle
global RBFtype; 
global RBFpar;
global RBFscale;
global A
% parameters for RBF 
soborder=2    % Sobolev space order for error estimation
RBFscale=1    % no artificial scaling
RBFtype='ms'  % Sobolev case 
RBFpar=soborder-1 % we are in dimension 2 
% points for interpolation
pts=[0 0 ; 1 0 ; 0 1];
% evaluation grid
he=0.05;
[xe, ye]=meshgrid(0:he:1,0:he:1);
pe=[xe(:), ye(:)];
imsk=find(pe(:,1)+pe(:,2)>1);
% mask for plotting only on the triangle
msk=ones(size(xe(:)));
msk(imsk)=0;
% the Lagrange basis functions for piecewise linear interpolation
u(:,1)=1-pe(:,1)-pe(:,2);
u(:,2)=pe(:,1);
u(:,3)=pe(:,2);
% for i=1:3
%   figure
%   surf(xe, ye, reshape(msk.*u(:,i),size(xe))) 
% end
% We now evaluate the generalized power function   
p0=frbf(0,0);
em=kermat(pe,pts); % kernel matrix on evaluation points 
a=kermat(pts,pts); % kernel matrix on data points 
pf=p0-2*sum(em'.*u')'+sum((u*a)'.*u')'; % generalized power function
figure
subplot(2,2,1)
surf(xe,ye,reshape(msk.*pf,size(xe)))
axis([0 1 0 1 0 0.3])
shading interp
% title('Power Function of PL Interpolant') 
% we now calculate the optimal power function.
% First, we need the Lagrange basis 
uk=(a\em')';
% then the power function 
pfk=p0-2*sum(em'.*uk')'+sum((uk*a)'.*uk')';
subplot(2,2,2)
surf(xe,ye,reshape(msk.*pfk,size(xe)))
shading interp
axis([0 1 0 1 0 0.3])
% title('Power Function of Optimal Interpolant') 
% then the difference
subplot(2,2,3)
surf(xe,ye,reshape(msk.*(pf-pfk),size(xe)))
shading interp
% title('Difference of Power Functions') 
