%--- Interpolates a function over the unit sphere using RBF-QR and
%--- RBF-Direct. The best result in each case is displayed and the
%--- max-errors are displayed as functions of the shape parameter. 
clear
N = 300; % The number of node points
Ne=10;   % The number of evaluation points in each dimension
d=3;     % The number of dimensions

%--- Approximately N Halton points in the unit sphere.
%--- Volume of a cube -1 to 1 is 2^3. Volume of sphere is pi r^3=pi.
xk_c=2*(haltonseq(round(N/pi*8),d)-0.5);
%
% Make into spherical coordinates with latitude longitude coordinates
% Order of spherical coordinates (r,th,phi)
%
[xk_p(:,3),xk_p(:,2),xk_p(:,1)]=cart2sph(xk_c(:,1),xk_c(:,2),xk_c(:,3));
xk_p(:,2) = pi/2-xk_p(:,2);

pos=find(xk_p(:,1)<=1); % Radius <= 1 in spherical coordinates
xk_c=xk_c(pos,:);  xk_p=xk_p(pos,:);

%--- Evaluation points from cartesian grid
t=-1:2/(Ne-1):1;
[x,y,z]=ndgrid(t,t,t);
x=x(:); y=y(:); z=z(:); 
pos = find(x.^2+y.^2+z.^2<=1);
x=x(pos);y=y(pos);z=z(pos);
[phi,th,r]=cart2sph(x,y,z);
th=pi/2-th;
xe_p=[r,th,phi];

%--- Right hand side and exakt solution  
fun = @(x,y,z) sin(x.^2+2*y.^2-z.^2)-sin(2*x.^2+(y-0.5).^2+(z+0.1).^2);
f = fun(xk_c(:,1),xk_c(:,2),xk_c(:,3));  
u = fun(x,y,z);  
%--- Interpolate and evaluate
epvec=logspace(-2,log10(2),25);
[cx,xx]=meshgrid(xk_c(:,1));
[cy,yy]=meshgrid(xk_c(:,2));
[cz,zz]=meshgrid(xk_c(:,3));
r2 = (xx-cx).^2 + (yy-cy).^2 + (zz-cz).^2;
[cx,xx]=meshgrid(xk_c(:,1),x);
[cy,yy]=meshgrid(xk_c(:,2),y);
[cz,zz]=meshgrid(xk_c(:,3),z);
r2e = (xx-cx).^2 + (yy-cy).^2 + (zz-cz).^2;
for k=1:length(epvec)
  ep=epvec(k)
  rbfqr(:,k) = RBF_QR_3D(ep,xk_p,xe_p,f);
  A = exp(-ep^2*r2);
  B = exp(-ep^2*r2e);
  rbfdir(:,k) = B*(A\f);
  err_rbfqr(k)  = max(abs(rbfqr(:,k) -u));
  err_rbfdir(k) = max(abs(rbfdir(:,k)-u)); 
end
%--- Display solutions and errors
figure(1),clf
loglog(epvec,err_rbfqr,'b',epvec,err_rbfdir,'r--')
axis tight
xlabel('\epsilon')
ylabel('Maximum error')
legend('RBF-QR','RBF-Direct')

