function hh = area(varargin)
%AREA  Filled area plot.
%   AREA(X,Y) produces a stacked area plot suitable for showing the
%   contributions of various components to a whole.  
%
%   For vector X and Y, AREA(X,Y) is the same as PLOT(X,Y) except that
%   the area between 0 and Y is filled.  When Y is a matrix, AREA(X,Y)
%   plots the columns of Y as filled areas.  For each X, the net
%   result is the sum of cooresponding values from the columns of Y.
%   X must be monotonic. 
%
%   AREA(Y) uses the default value of X=1:SIZE(Y,1).
%
%   AREA(X,Y,LEVEL) or AREA(Y,LEVEL) specifies the base level
%   for the area plot to be at the value y=LEVEL.  The default
%   value is LEVEL=0.
%
%   AREA(...,'Prop1',VALUE1,'Prop2',VALUE2,...) sets the specified
%   properties of the underlying patch objects.
%
%   H = AREA(...) returns a vector of handles to patch objects.
%
%   See also PLOT, BAR.

%   Clay M. Thompson 2-8-94
%   Copyright 1984-2000 The MathWorks, Inc. 
%   $Revision: 1.17 $  $Date: 2000/06/02 04:30:50 $

ax = newplot;
next = lower(get(ax,'NextPlot'));
hold_state = ishold;
nargs = nargin;

% Search for the beginning of the prop,value pairs.
for i=1:length(varargin),
  if isstr(varargin{i}), nargs = i-1; break, end
end

if nargs<3, level = 0; end

% Make sure x and y are the same size
if nargs<1, 
  error('Not enough input arguments.');
elseif nargs==1, % area(y)
  [msg,x,y] = xychk(varargin{1},'plot');
elseif nargs==2 % area(x,y) or area(y,level)
  % area(y,level)
  if ~isequal(size(varargin{1}),size(varargin{2})) & ...
      length(varargin{2})==1, 
    [msg,x,y] = xychk(varargin{1},'plot');
    level = varargin{2};
  else
    [msg,x,y] = xychk(varargin{1:2},'plot');
  end
else % area(x,y,level)
  [msg,x,y] = xychk(varargin{1:2},'plot');
  level = varargin{3};
end
if ~isempty(msg), error(msg); end
if all(size(level))~=1, error('LEVEL must be a scalar.'); end

if min(size(y))==1, y = y(:); x = x(:); end
[m,n] = size(y);

if n>1,
  % Check for the same x spacing
  if all(all(abs(diff(x')')<eps)), 
    y = cumsum(y')'; % Use fast calculation
  else
    xi = sort(x(:));
    yi = zeros(length(xi),size(y,2));
    for i=1:n,
      yi(:,i) = interp1(x(:,i),y(:,i),xi);
    end
    d = find(isnan(yi(:,1))); 
    if ~isempty(d), yi(d,1) = level(ones(size(d))); end
    d = find(isnan(yi)); 
    if ~isempty(d), yi(d) = zeros(size(d)); end
    x = xi(:,ones(1,n));
    y = cumsum(yi')';
    [m,n] = size(y);
  end
  xx = [x(1,:);x;flipud(x)];
  yy = [level(ones(m,1)) y];
  yy = [yy(1,1:end-1);yy(:,2:end);flipud(yy(:,1:end-1))];
else
  xx = [x(1,:);x;flipud(x)];
  yy = [level;y;level(ones(m,1))];
end

h = []; cc = ones(size(xx,1),1);
for i=1:size(y,2),
  h = [h,patch('xdata',xx(:,i),'ydata',yy(:,i),'cdata',i*cc, ...
      'faceColor','flat','edgecolor',get(ax,'xcolor'), ...
      varargin{nargs+1:end})];
end

if ~hold_state, 
   view(2); set(ax,'NextPlot',next); set(ax,'box','on')
   set(ax,'xlim',[min(xx(:)) max(xx(:))],'clim',[1 max(n,2)])
end

if nargout>0, hh = h; end
