function [A,B,x,y,z]=mquad_demo()
% MQUAD_DEMO Produces the Hardy multiquadric interpolant to the
% MATLAB's peaks function. To use the demo, enter
%
%             mquad_demo; 
%
% on the command line and press RETURN. A screen with cross hairs will
% appear. Select as many points as you like using the left mouse
% button. When you have finished, click with the right mouse
% button.

clf
axis([-3 3 -3 3]);
hold on
x=[];y=[];n=0;but=1;
while but==1,
[xi,yi,but]=ginput(1);
plot(xi,yi,'ro')
n=n+1;
y(n,1)=yi;x(n,1)=xi; 
end
hold off
x_cen=x; y_cen=y;
ydata=peaks(x_cen,y_cen);
ydata=[ydata;0];
N=length(x_cen);
[XX,YY]=rbf_inputs(x_cen,y_cen);
A=(1+XX.^2+YY.^2).^(1/2);
B=ones(N+1,N+1);
B(1:N,1:N)=A;
B(N+1,N+1)=0;
coefs=pinv(B)*ydata;
x_val=-3:0.1:3;
y_val=-3:0.1:3;
M=length(x_val);
[x,y]=meshgrid(x_val,y_val);
[xx,yy]=local_rbf_inputs(x_cen,y_cen,x(:),y(:));
z=reshape(coefs(1:N)'*(1+xx.^2+yy.^2).^(1/2) + coefs(N+1),M,M);
subplot(2,1,1), peaks(x,y)
subplot(2,1,2), surf(x,y,z), axis([-3 3 -3 3 -8 8]),...
    title('Multiquadric Interpolant')

function [XX,YY]=local_rbf_inputs(x_cen,y_cen,x_v,y_v)

%Compute XX(k,j)=x_v(k)-x_cen(j)
[XV,XC_T]=meshgrid(x_v,x_cen);
XX=XV-XC_T;

%Compute YY(k,j)=y_v(k)-y_cen(j)
[YV,YC_T]=meshgrid(y_v,y_cen);
YY=YV-YC_T;



