function zi = loess1d(xd,dd,xi,sx,check)
% loess filter estimate at xi of data [xd,dd] 
% zi = loess1d(xd,dd,xi,sx,check)
% xd is the coordinate
% dd is the data
% xi is the target coordinate
% sx is the loess filter width in units of xd
%
% If a 5th input (check) is given (it's value is not relevant) then
% the test for unsafe extrapolation in the case of a one-sided
% distribution of data within distance sx of the estimation 
% point is DISABLED
%
% John Wilkin

% preallocate output
zi = NaN*ones(size(xi));

% hold the indicies of valid data points - this will be clipped so that we
% only determine a smoothed value for points that were valid unsmoothed

datalocations = 1:length(dd);

% check for NaNs in the data
nodata = find(isnan(dd));
if ~isempty(nodata)
  dd(nodata) = [];
  xd(nodata) = [];
  datalocations(nodata) = [];
end
if isempty(dd)
  bell;disp('There are no data');return
end

% loop over index of input coordinates
for j=1:length(xi)
  
  xj = xi(j);

  % centre coordinates on x
  x = xd - xj;
  
  % distance metric
  r = abs(x/sx);
  
  % only use data within r<=1
  Q = find(r<=1);
  r = r(Q);
  x = x(Q);
  d = dd(Q);

  % convert to vectors
  x = x(:); r = r(:); d = d(:);

  if (min([length(find(x>0)) length(find(x<0))]) > 3) | nargin>4
  
    % compute local weighted least squares quadratic estimate (loess)
    
    % form matrix A of data coordinates
    A = [ones(size(x)) x x.^2];
    
    % calculate weights w
    w = (1-r.^3).^3;
    
    % apply weights
    A = w(:,ones([1 3])).*A;
    rhs = w.*d;
    
    % c is the vector of coefficients of the quadratic fit:
    % fit = c(1) + c(2)*x + c(3)*x^2
    
    % solve least-squares fit of A*c = d
    c = A\rhs;
    
    % evaluate fitted value:
    % formally, zj = [1 xj xj^2]*c, but since we moved
    % origin to xj the answer is just c(1)
    zj = c(1);
    
    zi(j) = zj;

  end
  
end