Wna = 1/4*(-2*yy(:,i2).*yy(:,i3)+yy(:,i2).^2+yy(:,i3).^2+xx(:,i2).^2 - ...
             2*xx(:,i2).*xx(:,i3)+xx(:,i3).^2).*zz; 
Wna(:) = Wna-1/16.*(yy(:,i2).^2-2*yy(:,i2).*yy(:,i3)+yy(:,i3).^2+xx(:,i2).^2- ...
         2.*xx(:,i2).*xx(:,i3)+xx(:,i3).^2).*(-xx(:,i2)+2.*xx(:,i1)-xx(:,i3)).*gx;
Wna(:) = Wna-1/16.*(yy(:,i2).^2-2.*yy(:,i2).*yy(:,i3)+yy(:,i3).^2+xx(:,i2).^2- ...
         2.*xx(:,i2).*xx(:,i3)+xx(:,i3).^2).*(-yy(:,i2)+2.*yy(:,i1)-yy(:,i3)).*gy;
Wna(:) = Wna./Area./len(:,i1);

Wnb = 1/4*(yy(:,i1).^2+yy(:,i1).*yy(:,i3)-3.*yy(:,i2).*yy(:,i1)+ ...
        3.*yy(:,i2).*yy(:,i3)-2.*yy(:,i3).^2+xx(:,i1).^2+xx(:,i1).*xx(:,i3)- ...
        3.*xx(:,i2).*xx(:,i1)+3.*xx(:,i2).*xx(:,i3)-2.*xx(:,i3).^2).*zz;
Wnb(:) = Wnb-1/16*(6*yy(:,i1).*xx(:,i2).*yy(:,i3)-3*yy(:,i1).^2.*xx(:,i2)- ...
         2*yy(:,i1).*xx(:,i1).*yy(:,i3)+2*yy(:,i1).^2.*xx(:,i1)- ...
         4*yy(:,i1).*xx(:,i3).*yy(:,i3)+yy(:,i1).^2.*xx(:,i3)- ...
         2*yy(:,i2).*xx(:,i3).*yy(:,i3)+2*yy(:,i2).*xx(:,i3).*yy(:,i1)+ ...
         2*yy(:,i2).*xx(:,i1).*yy(:,i3)-2*yy(:,i2).*xx(:,i1).*yy(:,i1)+ ...
         3*yy(:,i3).^2.*xx(:,i3)-3*yy(:,i3).^2.*xx(:,i2)-xx(:,i1).^2.*xx(:,i3)+ ...
         2*xx(:,i1).^3+10*xx(:,i2).*xx(:,i1).*xx(:,i3)-5*xx(:,i2).*xx(:,i1).^2- ...
         4*xx(:,i1).*xx(:,i3).^2-5*xx(:,i3).^2.*xx(:,i2)+3*xx(:,i3).^3).*gx;
Wnb(:) = Wnb-1/16*(-yy(:,i1).^2.*yy(:,i3)+2*yy(:,i1).^3+ ...
         10*yy(:,i2).*yy(:,i1).*yy(:,i3)-5*yy(:,i2).*yy(:,i1).^2- ...
         4*yy(:,i1).*yy(:,i3).^2-5*yy(:,i3).^2.*yy(:,i2)+3*yy(:,i3).^3+ ...
         6*yy(:,i2).*xx(:,i1).*xx(:,i3)-3*yy(:,i2).*xx(:,i1).^2- ...
         2*yy(:,i1).*xx(:,i1).*xx(:,i3)+2*yy(:,i1).*xx(:,i1).^2- ...
         4*yy(:,i3).*xx(:,i1).*xx(:,i3)+yy(:,i3).*xx(:,i1).^2- ...
         2*yy(:,i3).*xx(:,i2).*xx(:,i3)+2*yy(:,i3).*xx(:,i2).*xx(:,i1)+ ...
         2*yy(:,i1).*xx(:,i2).*xx(:,i3)-2*yy(:,i1).*xx(:,i2).*xx(:,i1)+ ...
         3*xx(:,i3).^2.*yy(:,i3)-3*yy(:,i2).*xx(:,i3).^2).*gy;
Wnb(:) = Wnb./Area./len(:,i2);

Wnc = 1/4*(yy(:,i2).*yy(:,i1)+yy(:,i1).^2-2*yy(:,i2).^2+3*yy(:,i2).*yy(:,i3)- ...
         3.*yy(:,i1).*yy(:,i3)+xx(:,i2).*xx(:,i1)+xx(:,i1).^2-2.*xx(:,i2).^2+ ...
         3.*xx(:,i2).*xx(:,i3)-3.*xx(:,i1).*xx(:,i3)).*zz;
Wnc(:) = Wnc-1/16*(yy(:,i1).^2.*xx(:,i2)-4*yy(:,i1).*xx(:,i2).*yy(:,i2)+ ...
         2*yy(:,i1).^2.*xx(:,i1)-2*yy(:,i2).*xx(:,i1).*yy(:,i1)- ...
         3*yy(:,i1).^2.*xx(:,i3)+6*yy(:,i2).*xx(:,i3).*yy(:,i1)- ...
         3*yy(:,i2).^2.*xx(:,i3)+3*yy(:,i2).^2.*xx(:,i2)+ ...
         2*yy(:,i1).*xx(:,i2).*yy(:,i3)-2*yy(:,i2).*xx(:,i2).*yy(:,i3)- ...
         2*yy(:,i1).*xx(:,i1).*yy(:,i3)+2*yy(:,i2).*xx(:,i1).*yy(:,i3)+ ...
         2*xx(:,i1).^3-xx(:,i2).*xx(:,i1).^2-4*xx(:,i2).^2.*xx(:,i1)- ...
         5*xx(:,i1).^2.*xx(:,i3)+10*xx(:,i2).*xx(:,i1).*xx(:,i3)+ ...
         3*xx(:,i2).^3-5*xx(:,i2).^2.*xx(:,i3)).*gx;
Wnc(:) = Wnc-1/16*(2*yy(:,i1).^3-yy(:,i2).*yy(:,i1).^2- ...
         4*yy(:,i2).^2.*yy(:,i1)-5*yy(:,i1).^2.*yy(:,i3)+ ...
         10*yy(:,i2).*yy(:,i1).*yy(:,i3)+3*yy(:,i2).^3- ...
         5*yy(:,i2).^2.*yy(:,i3)+yy(:,i2).*xx(:,i1).^2- ...
         4*yy(:,i2).*xx(:,i1).*xx(:,i2)+2*yy(:,i1).*xx(:,i1).^2- ...
         2*yy(:,i1).*xx(:,i2).*xx(:,i1)-3*yy(:,i3).*xx(:,i1).^2+ ...
         6*yy(:,i3).*xx(:,i2).*xx(:,i1)-3*yy(:,i3).*xx(:,i2).^2+ ...
         3*yy(:,i2).*xx(:,i2).^2+2*yy(:,i2).*xx(:,i1).*xx(:,i3)- ...
         2*yy(:,i2).*xx(:,i2).*xx(:,i3)-2*yy(:,i1).*xx(:,i1).*xx(:,i3)+ ...
         2*yy(:,i1).*xx(:,i2).*xx(:,i3)).*gy;
Wnc(:) = Wnc./Area./len(:,i3);

Wn = Wna(:,[1 2 3]) + Wnb(:,[3 1 2]) + Wnc(:,[2 3 1]);

% Find the nearest triangle (t)
t = tsearch(x,y,tri,xi,yi);

% Only keep the relevant triangles.
out = find(isnan(t));
if ~isempty(out), t(out) = ones(size(out)); end
tri = tri(t,:);
Area = Area(t,:);
len = len(t,:);
xx = xx(t,:);
yy = yy(t,:);
zz = zz(t,:);
gx = gx(t,:);
gy = gy(t,:);
gn = gn(t,:);
Wn = Wn(t,:);

% Compute Barycentric coordinates (w).  P. 78 in Watson.
w = 1/2.*((xx(:,i2)-repmat(xi,1,3)).*(yy(:,i3)-repmat(yi,1,3)) - ...
          (xx(:,i3)-repmat(xi,1,3)).*(yy(:,i2)-repmat(yi,1,3)))./Area;
w(out,:) = ones(length(out),3);

N1 = w(:,i1) + w(:,i1).^2.*w(:,i2) + w(:,i1).^2.*w(:,i3) - ...
               w(:,i1).*w(:,i2).^2 - w(:,i1).*w(:,i3).^2;
N2 = (xx(:,i2)-xx(:,i1)).*(w(:,i1).^2.*w(:,i2)+1/2.*w(:,i1).*w(:,i2).*w(:,i3))+ ...
     (xx(:,i3)-xx(:,i1)).*(w(:,i1).^2.*w(:,i3)+1/2.*w(:,i1).*w(:,i2).*w(:,i3));
N3 = (yy(:,i2)-yy(:,i1)).*(w(:,i1).^2.*w(:,i2)+1/2.*w(:,i1).*w(:,i2).*w(:,i3))+ ...
     (yy(:,i3)-yy(:,i1)).*(w(:,i1).^2.*w(:,i3)+1/2.*w(:,i1).*w(:,i2).*w(:,i3));
N1(out) = zeros(size(out));
N2(out) = zeros(size(out));
N3(out) = zeros(size(out));

M = 8*Area./len.*w(:,i1).*w(:,i2).^2.*w(:,i3).^2 ./ ...
                (w(:,i1)+w(:,i2)+(w(:,i1)+w(:,i2)==0)) ./ ...
                (w(:,i1)+w(:,i3)+(w(:,i1)+w(:,i3)==0));
M(out,:) = zeros(length(out),3);

zi = sum((N1.*zz + N2.*gx + N3.*gy + M.*(gn - Wn)).').';

zi = reshape(zi,siz);

if ~isempty(out), zi(out) = NaN; end
%------------------------------------------------------------

%------------------------------------------------------------
function zi = nearest(x,y,z,xi,yi)
%NEAREST Triangle-based nearest neightbor interpolation

%   Reference: David F. Watson, "Contouring: A guide
%   to the analysis and display of spacial data", Pergamon, 1994.

siz = size(xi);
xi = xi(:); yi = yi(:); % Treat these a columns
x = x(:); y = y(:); z = z(:); % Treat these as columns

% Triangularize the data
tri = delaunay(x,y,'sorted');
if isempty(tri), 
  warning('Data cannot be triangulated.');
  zi = repmat(NaN,size(xi));
  return
end

% Find the nearest vertex
k = dsearch(x,y,tri,xi,yi);

zi = k;
d = find(isfinite(k));
zi(d) = z(k(d));
zi = reshape(zi,siz);
%----------------------------------------------------------


%----------------------------------------------------------
function [xi,yi,zi] = gdatav4(x,y,z,xi,yi)
%GDATAV4 MATLAB 4 GRIDDATA interpolation

%   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.

xy = x(:) + y(:)*sqrt(-1);

% Determine distances between points
d = xy(:,ones(1,length(xy)));
d = abs(d - d.');
n = size(d,1);
% Replace zeros along diagonal with ones (so these don't show up in the
% find below or in the Green's function calculation).
d(1:n+1:prod(size(d))) = ones(1,n);

non = find(d == 0);
if ~isempty(non),
  % If we've made it to here, then some points aren't distinct.  Remove
  % the non-distinct points by averaging.
  [r,c] = find(d == 0);
  k = find(r < c);
  r = r(k); c = c(k); % Extract unique (row,col) pairs
  v = (z(r) + z(c))/2; % Average non-distinct pairs
  
  rep = find(diff(c)==0);
  if ~isempty(rep), % More than two points need to be averaged.
    runs = find(diff(diff(c)==0)==1)+1;
    for i=1:length(runs),
      k = find(c==c(runs(i))); % All the points in a run
      v(runs(i)) = mean(z([r(k);c(runs(i))])); % Average (again)
    end
  end
  z(r) = v;
  if ~isempty(rep),
    z(r(runs)) = v(runs); % Make sure average is in the dataset
  end

  % Now remove the extra points.
  x(c) = [];
  y(c) = [];
  z(c) = [];
  xy(c,:) = [];
  xy(:,c) = [];
  d(c,:) = [];
  d(:,c) = [];
  
  % Determine the non distinct points
  ndp = sort([r;c]);
  ndp(find(ndp(1:length(ndp)-1)==ndp(2:length(ndp)))) = [];

  warning(sprintf(['Averaged %d non-distinct points.\n' ...
       '         Indices are: %s.'],length(ndp),num2str(ndp')))
end

% Determine weights for interpolation
g = (d.^2) .* (log(d)-1);   % Green's function.
% Fixup value of Green's function along diagonal
g(1:size(d,1)+1:prod(size(d))) = zeros(size(d,1),1);
weights = g \ z(:);

[m,n] = size(xi);
zi = zeros(size(xi));
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
%----------------------------------------------------------