function [yi, varargout] = interp1_quick(x, y, xi, varargin)
% INTERP1_QUICK: a quick 1-dimensional interpolation - useful for taking slices
%
%    INPUT:
% X: an MxN  matrix of values; each column must be strictly monotonic
% Y: an MxN matrix
% XI: either a scalar or a vector of length N
% VARARGIN: further MxN matrices that will be used for interpolation in the
%           same way that Y is.
%
%    OUTPUT:
% YI: a vector of length N. For the I'th column in X the position of X(I) is
%     found by linear interpolation. YI is the corresponding value in Y.
% VARARGOUT: Further output vectors of length N. In the same way that YI
%            results from interpolation in Y these result from interpolation
%            in the matrices in VARGIN
%
%   NOTES
%
% 1) The same value would be returned by:
%     yii = zeros(1, n_cols);
%     n_cols = size(x, 2);
%     for ii = 1:n_cols
%       yii(ii) = interp1(x(:, ii), y(:, ii), xi(ii));
%     end
% However this is version is significantly faster - by a factor of about
% 10 for N ~ 50 and 100 for N ~ 500.
%
% 2) If the input matrix X is not strictly monotonic in each column then
% empty matrices are returned for every output argument as well as an error
% message being displayed
%
% Tested by running t_interp1_quick.m

% $Id: interp1_quick.m,v 1.3 1998/05/11 01:58:04 mansbrid Exp $
% Copyright J. V. Mansbridge, CSIRO, Thu May  7 12:47:30 EST 1998

si_x = size(x);
n_rows = si_x(1);
n_cols = si_x(2);
si_y = size(y);
si_xi = size(xi);

% Check that x and y are both matrices of the same dimensions.

if length(si_y) ~= 2
  error('Y must be a matrix')
elseif sum(abs(si_y - si_x)) ~= 0
  error('X and Y must be the same size')
end

% Determine type of montonicity and check for non-monotonic behaviour.

dd = x(2:end, :) - x(1:end-1, :);
min_dd = min(dd(:));
max_dd = max(dd(:));
if min_dd > 0
  mono_inc = 1;
elseif max_dd < 0
  mono_inc = 0;
else
  disp('interp1_quick: The columns of X are not strictly monotonic')
  yi = [];
  if nargout > 1
    for ii = 1:nargout - 1
      varargout{ii} = [];
    end
  end
  return
end

% If xi is a scalar then use this for value to be interpolated to. If it is a
% vector then the comparison with x must be with a matrix or vector according
% to the context.

if max(si_xi) == 1
  xi_matrix = xi;
  xi_is_scalar = 1;
else
  xi_vector = xi(:);
  xi_matrix = repmat(xi_vector, 1, n_rows)';
  xi_vector = xi_vector';
  xi_is_scalar = 0;
end  

% For each column find the value index such that xi(ii) is bounded by
% x(index(ii), ii) and x(index(ii)+1, ii). It is assumed that each x(:, ii)
% column is monotonic of the same type.

if mono_inc
  index = sum(x < xi_matrix, 1);
else
  index = sum(x > xi_matrix, 1);
end

% Find the unacceptable indices, record them, and then replace them with
% acceptable values.

ff_low = find(index == 0);
len_low = length(ff_low);
if (len_low > 0)
  index(ff_low) = repmat(1, size(ff_low));
end
ff_up = find(index == n_rows);
len_up = length(ff_up);
if (len_up > 0)
  index(ff_up) = repmat((n_rows - 1), size(ff_up));
end

% Find the interpolated values yi.

index_full = index + n_rows*(0:(n_cols - 1));
index_full_p = index_full + 1;
x_diff = x(index_full_p) - x(index_full);
if xi_is_scalar
  diff_1 = x(index_full_p) - xi;
  diff_2 = xi - x(index_full);
else
  diff_1 = x(index_full_p) - xi_vector;
  diff_2 = xi_vector - x(index_full);
end
yi = (y(index_full).*diff_1 + y(index_full_p).*diff_2)./x_diff;

% replace the bad values (extrapolations) with NaNs.

if (len_low > 0)
  yi(ff_low) = repmat(NaN, size(ff_low));
end
if (len_up > 0)
  yi(ff_up) = repmat(NaN, size(ff_up));
end

% Now deal with the extra input and output variables.

num_extra = min([(nargout-1) (nargin -3)]);
if num_extra > 0
  for ii = 1:num_extra

    % Check that x and y_tmp are both matrices of the same dimensions.

    y_tmp = varargin{ii};
    si_y_tmp = size(y_tmp);
    if length(si_y_tmp) ~= 2
      error('Y_TMP must be a matrix')
    elseif sum(abs(si_y_tmp - si_x)) ~= 0
      error('X and Y_TMP must be the same size')
    end

    % Find the interpolated values yi_tmp.

    yi_tmp = (y_tmp(index_full).*diff_1 + y_tmp(index_full_p).*diff_2)./x_diff;

    % replace the bad values (extrapolations) with NaNs.

    if (len_low > 0)
      yi_tmp(ff_low) = repmat(NaN, size(ff_low));
    end
    if (len_up > 0)
      yi_tmp(ff_up) = repmat(NaN, size(ff_up));
    end
    varargout{ii} = yi_tmp;
  end
end