plotstack(X,varargin)

function [ax,ls] = plotstack(X,varargin)
% PLOTSTACK plots the columns in a matrix, one above the other.
%
% CALL
%   [ax,ls] = plotstack(X,...);           % X is the 2D array 
%   [ax,ls] = plotstack(X,Y,...);         % Y is the 2D array, X is the x-axis 
%   [ax,ls] = plotstack(X,plotfun,...);   % X is the 2D array, plotfun is
%                                   % plotting function
%   [ax,ls]ls = plotstack(X,Y,plotfun,...);
% 
% INPUT
%   plotfun: (function handle) plotting function to use. E.g. @plot, @stem,
%            @stairs. Default is @plot.
%   ...:     plotfun-specific properties and values, such as 'MarkerSize', 
%            'Linewidth' etc.
%            An additional property is 'Vertical' to plot from
%            to top to bottom (see examples)
% 
% OUTPUT
%   ax: axis handles (for futher manipulation).
%   ls: lineseries handles (for further manipulation).
%
% EXAMPLES
%   x = randn(100,10);  % make random data
%   figure, plotstack(x,'k');
%   figure, plotstack(1:100,x,'m-*');
%   figure, plotstack(x,@stem,'MarkerSize',2); 
%   figure, ax = plotstack((1:50)/50,x(1:50,1:5), @stairs, 'r','Vertical', 'Linewidth',2);
%           xlabel('time'), 
%           for i = 1:length(ax), 
%               set(get(ax(i),'Ylabel'),'String',sprintf('signal %d',i)); 
%           end
%
% C. Saragiotis, June 2011


%% input check
if nargin == 1
    t = 1:size(X,1);
    plotfun = @plot;
elseif nargin == 2
    if isnumeric(varargin{1})
        t = X;
        X = varargin{1};
        plotfun = @plot;
        varargin(1) = [];
    elseif isa(varargin{1},'function_handle')
        t = 1:size(X,1);
        plotfun = varargin{1};
        varargin(1) = [];
    elseif ischar(varargin{1})
        t = 1:size(X,1);
        plotfun = @plot;
    end
elseif nargin>2
   if isnumeric(varargin{1})
       t = X;
       X = varargin{1};
       if isa(varargin{2},'function_handle')
           plotfun = varargin{2};
           varargin(1:2) = [];
       else
           plotfun = @plot;
           varargin(1) = [];
       end
   elseif isa(varargin{1},'function_handle')
       t = 1:size(X,1);
       plotfun = varargin{1};
       varargin(1) = [];
   else
       t = 1:size(X,1);
   end
end

if isvector(X), error('%s: X is scalar'); end

[nt,nx] = size(X);

if length(t)~=nt
   error('%s: dimensions of X and Y don''t match.',upper(mfilename));
end

nvarg  = length(varargin);
isvert = 0;
if nvarg > 0
    for i = 1:nvarg
        if strcmpi(varargin{i},'Vertical')
            isvert = 1;
            varargin(i) = [];
            break;
        end
    end    
end

ax = zeros(nx,1);
ls = zeros(nx,1);

%% Main
switch isvert
    case 0
        for ix = 1:nx
            ax(ix) = subplot(nx,1,ix);
            ls(ix) = plotfun(t,X(:,ix),varargin{:});

            set(gca,'YTickLabel','');
            set(gca,'YTick',[]);

            set(gca,'Box','off');        
            if ix<nx
                set(gca,'XTickLabel','');
                set(gca,'XTick',[]);
            end
            if ix == 1,
                pos = get(gca,'Position');
                top = pos(2) + pos(4);
            elseif ix == nx
                pos = get(gca,'Position');
                bot = pos(2);
                hgt = (top-bot)/nx;
            end
        end 
        
        for ix = 1:nx
            pos = get(ax(ix),'Position');
            pos(2) = top - ix*hgt ;
            pos(4) = hgt;
            set(ax(ix),'ActivePositionProperty','Position')
            set(ax(ix),'Position',pos)
        end
        
    case 1
        for ix = nx:-1:1 % inverse count, so that CurrentAxes 
                         % is the leftmost, to set e.g. the xlabel
            ax(ix) = subplot(1,nx,ix);
            ls(ix) = plotfun(t,X(:,ix),varargin{:});
            view([90,90])

            set(gca,'YTickLabel','');
            set(gca,'YTick',[]);

            set(gca,'Box','off');        
            if ix>1
                set(gca,'XTickLabel','');
                set(gca,'XTick',[]);
            end
            if ix == nx
                pos = get(gca,'Position');
                rgt = pos(1)+pos(3);
            elseif ix == 1,
                pos = get(gca,'Position');
                lft = pos(1);
                wdt = (rgt-lft)/nx;
            end
        end
   
        for ix = nx:-1:1
            pos = get(ax(ix),'Position');
            pos(1) = lft + (ix-1)*wdt;
            pos(3) = wdt;
            set(ax(ix),'ActivePositionProperty','Position')
            set(ax(ix),'Position',pos)
        end

end