
function [ixmin, ixmax, iymin, iymax] = inside1(x, y, ax)

% Given a plot of array x vs. vector y, and a coordinate pair
% defining a rectangular subset of the plot
% ax = [xmin xmax ymin ymax], it returns a vector iax =
% [ixmin ixmax iymin iymax] where 
% (ixmin,iymin) is the point nearest and >= (ax(1),ax(3)) 
% (ixmax,iymax) is the point nearest and <= (ax(2),ax(4))
% Hence, it does not necessarily include all the points that
% fall in the rectangle, nor does it exclude points that
% fall outside the rectangle--it all depends on the shape of the profiles.

% set to default values:
ixmin = [];
ixmax = [];
iymin = [];
iymax = [];

irow = find(y >= ax(3) & y <= ax(4));
if ~isempty(irow)
  iymin = min(irow);  % assumes y is already sorted 
  iymax = max(irow);
  ixmin = find(x(iymin,:) >= ax(1));
  if ~isempty(ixmin)
    ixmin = ixmin(1);
    ixmax = find(x(iymax,:) <= ax(2));
    if ~isempty(ixmax)
      ixmax = ixmax(length(ixmax));
    end
  end
end
