function [xi,yi,zi] = griddat(x,y,z,xi,yi,fun,n)
% GRIDDAT Data gridding.
%	ZI = GRIDDAT(X,Y,Z,XI,YI) returns matrix ZI containing elements
%	corresponding to the elements of matrices XI and YI and determined by
%	interpolation within the 2-D function described by the (usually)
%	nonuniformly-spaced vectors (X,Y,Z).
%	XI can be a row vector, in which case it specifies a matrix with
%	constant columns. Similarly, YI can be a column vector and it 
%	specifies a matrix with constant rows.
%
%	[XI,YI,ZI] = GRIDDAT(X,Y,Z,XI,YI) returns the XI and YI formed
%	this way, which are the same as the matrices returned by MESHGRID.
%
%	[...] = GRIDDAT(X,Y,Z,XI,YI,FUN,N)
%	uses function name or expression FUN to calculate Green's
%	functions usen in interpolation. It should both accept
%	and return distance variable 'r', like FUN='exp(-r^2)'.
%	Default setting is equivalent to FUN='r^2*log(r)'.
%	N is maximum number of points allowed in an elementary
%	block in QUADTREE division procedure.
%
%	GRIDDAT uses a variation of objective mapping method.
%	It divides the region of interest in smaller regions
%	using QUADTREE and for each interpolation point uses
%	only the points in the same and neighbouring regions.

%  Kirill K. Pankratov,  kirill@plume.mit.edu
%  01/31/95

n_dflt = 25;    % Default for max. number of points
                % in elementary block

 % Handle input .......................
error(nargchk(5,7,nargin))
[msg,x,y,z,xi,yi] = xyzchk(x,y,z,xi,yi);
if length(msg)>0, error(msg); end

if nargin<7, n = n_dflt; end
is_call = 0;
if nargin==6
  if isstr(fun)
    is_call=1;
  else
    n = fun(1);
  end
end
n = max(n,1);

 % Set up function call .....................................
if is_call
  is_fun = exist(fun); % Check if a function or expression
  is_fun = is_fun & any(isletter(fun));
  if is_fun          % .. Function
    call = ['r=' fun '(r)'];

  else               % .. Expression

    % Change to array operations (.* , .^, ./) :
    n_ins = fun=='*' | fun=='/' | fun=='^';
    nn = find(n_ins);
    n_ins = cumsum(n_ins+1);
    call = zeros(size(max(n_ins)));
    call(n_ins) = fun;
    call(~call) = '.'*ones(size(nn));
    % Add output and semicolon
    call = setstr(['r=' call ';']);
  end
end

szi = size(xi);
l  = prod(size(x));
li = prod(szi);
if l<=100, n=100; end

 % Combine known and interp. points in one vector
x = [x(:); xi(:)];
y = [y(:); yi(:)];
z = [z(:); zeros(li,1)];
s = [ones(l,1); zeros(li,1)];


 % Calculate "quadtree"
[ind,bx,by,Nb] = quadtree(x,y,s,n);
if Nb==[], Nb=1; end
lb = size(Nb,1);

 % Construct sparse matrix for indexing
I = sparse(ind,1:length(ind),ind);

 % For each block and its neighbours calculate
 % inverse Green's function matrix and interpolation
for jj = 1:lb

  % Neighbours of a current block
  i_nb = find(Nb(jj,:));

  % Current points
  if length(i_nb)<=1
    nn = full(find(I(i_nb,:)));
  else
    nn = full(find(any(I(i_nb,:))));
  end      % All pts. in cur. blocks

  nni = full(find(I(jj,:)));
  nni = nni(~s(nni)); % Interp. points
  nn = nn(s(nn));     % Basis points

  ll  = length(nn);
  lli = length(nni);
  ol  = ones(1,ll);
  oli = ones(1,lli);

  xx = x(nn); yy = y(nn);
  D = xx(:,ol);
  r = (D-D').^2;
  D = yy(:,ol);
  r = r+(D-D').^2;
  mask = find(~r);

  if is_call
    r = sqrt(r);
    eval(call);
  else
    r(mask) = ones(size(mask));
    r = r.*log(r);
  end

  is_1 = size(r,1)<=1;
  if ~is_1
    w = r \ z(nn);  % Weights
  end

  % Interpolated points
  D = xx(:,oli)'-x(nni,ol);
  r = D.^2;
  D = yy(:,oli)'-y(nni,ol);
  r = r+D.^2;
  mask = find(~r);

  if is_call
    r = sqrt(r);
    eval(call);
  else
    if mask~=[]
      r(mask) = ones(size(mask));
    end
    r = r.*log(r);
  end

  if ~is_1
    z(nni) = r*w;
  end

end

zi = z(l+1:l+li);    % Extract interpolated values
zi = reshape(zi,szi(1),szi(2));

if nargout<=1,
  xi = zi;
end
