function [r1,r2] = bins(A, edges, areaflag)
% BINS  List count of items in bins
%   C = BINS(A,EDGES) Returns a list of counts, indicating how many
%   elements in the array A are in each bin, whose boundaries are given
%   in EDGES.  To be precise, the number of elements in bin i is the number
%   of values in A satisfying  EDGES(i) < A <= EDGES(i+1).
%   Use -inf and/or inf as the first or last edge if you want
%   the first or last bin to be unbounded.
%   Note that the number of bins is one less than the number of edges.
%   C = BINS(A,EDGES,1) rescales the counts, so that they can be correctly
%   plotted on a probability histogram.  That is, the counts are each divided
%   by the total number of elements in A, and then the count in each bin is
%   rescaled by dividing by the width of that bin, so that the counts are
%   in fact densities.
%   Finally, [E,C] = BINS(A,EDGES) returns the counts (in C) along with
%   the midpoints of the bins (in E); these are suitable for then doing
%   e.g. BAR(E,C) to make a histogram.  For bins which have +/- infinity as
%   one of their edges, the midpoint of that bin will be the same distance
%   from the (finite) edge of the bin as the midpoint of the next interval is.
%
%   If EDGES is a positive scalar N, the data will be divided into N bins,
%   where the range of each bin is of width (max-min)/N using the min and
%   max of the data.
%   If EDGES is a negative scalar -N, the data will be divided into N bins,
%   where each bin has roughly the same number of items.
%
%   Note: you can get some similar functionality from HIST and/or HISTC.
%
% Written by David Hiebeler, http://www.math.umaine.edu/faculty/hiebeler


if (nargin == 2)
  areaflag = 0;
end

tmpA = A(:);  % if A is multidimensional, convert it into a 1D column vector
retval = [];  % start with empty array

if (ndims(edges) ~= 2)
  error('EDGES must be a vector')
end
if ((size(edges,1) ~= 1) & (size(edges,2) ~= 1))
  error('EDGES must be a vector')
end

edgeslen = max(size(edges));

% divide into n equally-sized bins (i.e. the max val minus min val in each
% bin is the same)
if ((edgeslen == 1) & (edges > 0))
  numbins = edges;
  binsize = (max(tmpA) - min(tmpA))/numbins;
  epsilon = binsize/1000;
  startval = min(tmpA);
  edges = startval-epsilon;
  for i=1:numbins-1
    edges = [edges startval+i*binsize];
  end
  edges = [edges max(tmpA)];
  edgeslen = numbins+1;
end


% divide into n bins which contain roughly the same number of elements
if ((edgeslen == 1) & (edges < 0))
  numbins = -edges;
  sA = sort(tmpA);
  len = size(sA, 1);
  numperbin = len/edges;
  edges = -inf;
  for i=1:numbins-1
    tmpindex = floor(i*len/numbins);
    edges = [edges sA(tmpindex)];
  end
  edges = [edges inf];
  edgeslen = numbins+1;
end

for i=1:edgeslen-1
  retval(i) = sum((tmpA > edges(i)) & (tmpA <= edges(i+1)));
end

if (areaflag)
  retval = retval ./ diff(edges) / length(tmpA);
end

if (nargout == 1)
  r1 = retval;
elseif (nargout == 2)
  r2 = retval;
  r1 = (edges(2:end) + edges(1:end-1))/2;
  if (edges(1) == -inf)
    r1(1) = edges(2) - (r1(2)-edges(2));
  end
  if (edges(end) == inf)
    r1(end) = edges(end-1) + (edges(end-1)-r1(end-1));
  end
end