function [xi,yi,zi] = griddata(x,y,z,xi,yi)
%GRIDDATA Data gridding.
%	ZI = GRIDDATA(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] = GRIDDATA(X,Y,Z,XI,YI) returns the XI and YI formed
%	this way, which are the same as the matrices returned by MESHGRID.
%
%	GRIDDATA uses an inverse distance method.
%
%	See also INTERP2, INTERP1.

%	Copyright (c) 1984-1998 by The MathWorks, Inc.

%       Reference:  David T. Sandwell, Biharmonic spline
%       interpolation of GEOS-3 and SEASAT altimeter
%       data, Geophysical Research Letters, 2, 139-142,
%       1987.  Describes interpolation using value or
%       gradient of value in any dimension.

% MATLAB 5 makes use of 'method', ignore it for now.
error(nargchk(5,5,nargin)) 

[msg,x,y,z,xi,yi] = xyzchk(x,y,z,xi,yi);

if length(msg)>0, error(msg); end
xy = x(:) + y(:)*sqrt(-1);

% Determine weights for interpolation
d = xy * ones(1,length(xy));
d = abs(d - d.');

mask = find(d == 0);
if length(mask) ~= length(d),
  error('All the (X,Y) points must be distinct.');
end

d(mask) = ones(length(mask),1);
g = (d.^2) .* (log(d)-1);   % Green's function.
g(mask) = zeros(length(mask),1); % Value of Green's function at zero
weights = g \ z(:);

[m,n] = size(xi);
zi = zeros(m,n);
jay = sqrt(-1);
xy = xy.';

% Evaluate at requested points (xi,yi).  Loop to save memory.
for i=1:m
  for j=1:n
    d = abs(xi(i,j)+jay*yi(i,j) - xy);
    mask = find(d == 0);
    if length(mask)>0, d(mask) = ones(length(mask),1); end
    g = (d.^2) .* (log(d)-1);   % Green's function.
    % Value of Green's function at zero
    if length(mask)>0, g(mask) = zeros(length(mask),1); end
    zi(i,j) = g * weights;
  end
end

if nargout<=1,
  xi = zi;
end

function [msg,x,y,z,out5,out6] = xyzchk(arg1,arg2,arg3,arg4,arg5)
%XYZCHK Check arguments to 3-D data routines.
%   [MSG,X,Y,Z,C] = XYZCHK(Z), or
%   [MSG,X,Y,Z,C] = XYZCHK(Z,C), or
%   [MSG,X,Y,Z,C] = XYZCHK(X,Y,Z), or
%   [MSG,X,Y,Z,C] = XYZCHK(X,Y,Z,C), or
%   [MSG,X,Y,Z,XI,YI] = XYZCHK(X,Y,Z,XI,YI) checks the input aguments
%   and returns either an error message in MSG or valid X,Y,Z (and
%   XI,YI) data.

%   Copyright (c) 1984-1998 by The MathWorks, Inc.
%   $Revision: 1.4 $  $Date: 1997/12/12 20:18:25 $

error(nargchk(1,6,nargin));

msg = [];
out5 = []; out6 = [];
if nargin==1, % xyzchk(z)
  z = arg1;
  if isstr(z),
    msg = 'Input arguments must be numeric.';
    return
  end
  [m,n] = size(z);
  [x,y] = meshgrid(1:n,1:m);
  out5 = z; % Default color matrix
  return

elseif nargin==2, % xyzchk(z,c)
  z = arg1; c = arg2;
  [m,n] = size(z);
  [x,y] = meshgrid(1:n,1:m);
  if ~isequal(size(z),size(c)),
    msg = 'Matrix C must be the same size as Z.';
    return
  end
  out5 = c;
  return

elseif nargin>=3, % xyzchk(x,y,z,...)
  x = arg1; y = arg2; z = arg3;
  [m,n] = size(z);
  if ~isvector(z), % z is a matrix
    % Convert x,y to row and column matrices if necessary.
    if isvector(x) & isvector(y),
      [x,y] = meshgrid(x,y);
      if size(x,2)~=n & size(y,1)~=m,
        msg = 'The lengths of X and Y must match the size of Z.';
        return
      elseif size(x,2)~=n,
        msg = 'The length of X must match the number of columns of Z.';
        return
      elseif size(y,1)~=m,
        msg = 'The length of Y must match the number of rows of Z.';
        return
      end
    elseif isvector(x) | isvector(y),
      msg = 'X and Y must both be vectors or both be matrices.';
      return
    else
      if ~isequal(size(x),size(y)) | ~isequal(size(x),size(z)),
        msg = 'Matrices X and Y must be the same size as Z.';
        return
      end
    end
  else % z is a vector
    if ~isvector(x) | ~isvector(y),
      msg = 'When Z is a vector, X and Y must also be vectors.';
      return
    elseif (length(x)~=length(z) | length(y)~=length(z)) & ...
           ~((length(x)==size(z,2)) & (length(y)==size(z,1)))
      msg = sprintf(['X and Y must be same length as Z or the lengths \n', ...
                    'of X and Y must match the size of Z.']);
      return
    end
  end
end

if nargin==4, % xyzchk(x,y,z,c)
  c = arg4;
  if ~isequal(size(z),size(c))
    msg = 'Matrix C must be the same size as Z.';
    return
  end
  out5 = c;
  return
end

if nargin==5, % xyzchk(x,y,z,xi,yi)
  xi = arg4; yi = arg5;

  if automesh(xi,yi),
    [xi,yi] = meshgrid(xi,yi);
  elseif ~isequal(size(xi),size(yi)),
    msg = 'XI and YI must be the same size or vectors of different orientations.';
  end
  out5 = xi; out6 = yi;
end

function tf = isvector(x)
%ISVECTOR True if x has only one non-singleton dimension.
tf = (length(x) == prod(size(x)));
