% Inverse distance squared interpolation
% v(i,j,k) = (f(1)/ds(1) + f(2)/ds(2) + ...)/(1/ds(1) + 1/ds(2) + ...)
% where f(p) is the value of the pth data point, ds(p) is the square of 
% the distance between the pth data point and the i,j,k grid point plus 
% a small distance epsilon
function v = inv_dist_squ_interp(x,y,z,dpts);
   [npts ncols] = size(dpts);
   v = zeros(length(y), length(x), length(z));

   % Compute v(i,j,k) for all i,j,k
   for i=1:length(y)              % x coordinate
      for j=1:length(x)           % y coordinate
         for k=1:length(z)        % z coordinate
            dx2 = (x(j)-dpts(:,1)).^2;
            dy2 = (y(i)-dpts(:,2)).^2;
            dz2 = (z(k)-dpts(:,3)).^2;
            dp = 1./(eps + dx2 + dy2 + dz2);
            fp = dpts(:,4)./(eps + dx2 + dy2 + dz2);
            v(i,j,k) = sum(fp)/sum(dp);
         end
      end
   end
end