% interp_test.m: Matlab program to test interp_shell.c, a Matlab
%                mex file


%mex interp_shell.c

% set up the grids where the function to be
% interpolated, v.m, is to be evaluated

xgrid=(0:.1:1);
nx=size(xgrid,2);

ygrid=(0:.1:1);
ny=size(ygrid,2);

zgrid=(0:.1:1);
nz=size(zgrid,2);

grid=[xgrid ygrid zgrid]';
n=nx*ny*nz;
vf=zeros(n,1);
p=zeros(n,1);

% now do a 3 level do loop to evaluate the function
% at those grid points

l=1;

for i=1:nx;
  for j=1:ny;
    for k=1:nz;

    vf(l)=v(xgrid(i),ygrid(j),zgrid(k));

    l=l+1;

    end;
  end;
end;

for i=1:nx;
   grid(i)=xgrid(i);
end;

for j=1:ny;
   grid(j+nx)=ygrid(j);
end;

for k=1:nz;
    grid(k+nx+ny)=zgrid(k);
end;

ns=1000;
s=10*rand(ns,3);
iv=zeros(ns,1);
tv=zeros(ns,1);

for i=1:ns;
  tv(i)=v(s(i,1),s(i,2),s(i,3));
end;

vs=s';
vs=vs(:);

gridnum=[nx ny nz]';

intmeth=0;
interp_shell(ns,3,grid,gridnum,vs,vf,p,iv,intmeth);

if (intmeth);
   fprintf('Running multidimensional interpolation routine, simplicial interpolation\n');
else;
   fprintf('Running multidimensional interpolation routine, multilinear interpolation\n');
end;
	

fprintf('max difference between interpolated and true v: %9.6e\n',max(abs(iv-tv)));
