function [zi,S,A,W]=scattered_interp2(x,y,z,e,xi,yi,extrap)
%scattered_interp2::2D linear interpolation from ungridded data with prespecified triangulation. 
%function zi=scattered_interp2(x,y,z,e,xi,yi)
%Hacked from griddata to allow user specification of the triangulation
%LINEAR Triangle-based linear interpolation or
%
%Inputs:
%       x(Nx1) x coordinates of data points
%       y(Nx1) y coordinates of data points
%       z(Nx1) z values of data.
%       e(Ntrix1) planar triangulation of x,y data points(such as e=delaunay(x,y) )
%       xi(Npx1) x coordinates of points to interpolate data to.
%       yi(Npx1) y coordinates of points to interpolate data to.
%       extrap: 1 to extrapolate to data outside the triangulation 0 not
%       to.
%
%Outputs:
%       zi(Npx1) estimate of z at xi yi.
%
%   Reference: David F. Watson, "Contouring: A guide
%   to the analysis and display of spacial data", Pergamon, 1994.
%
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%
% Keston Smith & Ata Bilgili May 2001
% Update by Ata Bilgili, Dartmouth College, 12 December 2002.
%   - Added warning about Matlab's delaunay.m routine which gets into effect when e is empty.
%
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
if nargin<7,
    extrap=1
end
siz = size(xi);
xi = xi(:); yi = yi(:); % Treat these as columns
x = x(:); y = y(:); % Treat these as columns

% Triangularize the data using Matlab's delaunay algorithm.

if isempty(e),
    e = delaunay(x,y);
    disp('The triangulation array e is empty!');
    disp('Using Matlabs delaunay.m routine for triangulation in scattered_interp2.m')
    disp('Make sure this triangulation is OK if you encounter bathymetry anomalies in the final grid')
end

if isempty(e),
  warning('Data cannot be triangulated.');
  zi = repmat(NaN,size(xi));
  return
end
% Find the nearest triangle (t)
t = tsearchsafe(x,y,e,xi,yi);
[n,m]=size(z);
zi=nan*ones(size(xi),m);
for k=1:length(t),
    if ~isnan(t(k))
        j=t(k);
        X=[ones(3,1),x(e(j,:)'),y(e(j,:)')];
        Z=z(e(j,:)',:);
        a=X\Z;
        zi(k,:)=a(1,:)+a(2,:)*xi(k)+a(3,:)*yi(k);
    else
        if extrap==1,
            disp('points outside of data tringulation found![or error in tsrchmx] extrapolating from nearest triangle');
            [m,j]=min(abs(mean(   x(e')-xi(k) + i*(y(e')-yi(k))   ))); 
        X=[ones(3,1),x(e(j,:)'),y(e(j,:)')];
        Z=z(e(j,:)');
        a=X\Z;
        zi(k,:)=a(1,:)+a(2,:)*xi(k)+a(3,:)*yi(k);
        end
    end        
end
