function msubplot(varargin)
% MSUBPLOT Replacement for Matlab's subplot funcion
%
% Usage:
% 1. Set up the subplot grid
%       msubplot(m, n, ...);
%
% 2. Create plots
%       msubplot(1);
%       plot(...);
%       msubplot(2);
%       plot(...);
%    etc.
%
% Inputs:
%    m, n:        Number of plots in grid
%
% Optional arguments
%    Border,                   : Border around subplots
%    BorderLeft, BorderRight   : Border left / right of subplots 
%    BorderTop, BorderBottom   : Border above / below subplots
%
%    Spacing,                  : Spacing between subplots
%    SpacingX, SpacingY        : Spacing in x or y direction btwn. subplots
%    
%    Heights, Widths           : Relative heights / widths of subplots
%                                (vector of length m or n respectively).
%
% Copyright 2016, Mikkel N. Schmidt, Technical University of Denmark

% Parse arguments
% -------------------------------------------------------------------------
% Input arguments
if nargin>=3 && isnumeric(varargin{1}) && isnumeric(varargin{2}) && isnumeric(varargin{3})
    opts = mgetopt(varargin(4:end));
    opts.m = varargin{1};
    opts.n = varargin{2};
    opts.p = varargin{3};
elseif nargin>=2 && isnumeric(varargin{1}) && isnumeric(varargin{2})
    opts = mgetopt(varargin(3:end));
    opts.m = varargin{1};
    opts.n = varargin{2};
elseif nargin>=1 && isnumeric(varargin{1})
    opts = mgetopt(varargin(2:end));
    opts.p = varargin{1};
else
    opts = mgetopt(varargin);
end

% Get options from UserData
uopts = get(gcf, 'UserData');
if isstruct(uopts)
    % Merge in UserData options
    FieldNames = fieldnames(uopts);
    for k = 1:length(FieldNames)
        % Don't overwrite new options
        if ~isfield(opts, FieldNames{k})
            opts.(FieldNames{k}) = uopts.(FieldNames{k});
        end
    end
end

% Parse optional arguments
p = mgetopt(opts, 'p', 1);
n = mgetopt(opts, 'n', 1);
m = mgetopt(opts, 'm', 1);
% Convert to linear index if plot number is given as a vector
if length(p)==2, p = p(2)+n*(p(1)-1); end
    
% Borders around plots
Border = mgetopt(opts, 'Border', 0.1);
BorderLeft = mgetopt(opts, 'BorderLeft', Border);
BorderRight = mgetopt(opts, 'BorderRight', Border);
BorderTop = mgetopt(opts, 'BorderTop', Border);
BorderBottom = mgetopt(opts, 'BorderBottom', Border);

% Spacing between subplots
Spacing = mgetopt(opts, 'Spacing', 0.1);
SpacingX = mgetopt(opts, 'SpacingX', Spacing);
SpacingY = mgetopt(opts, 'SpacingY', Spacing);

% Relative heights and widths
Heights = mgetopt(opts, 'Heights', ones(m, 1));
Widths = mgetopt(opts, 'Widths', ones(n, 1));

% Axis handles
H = mgetopt(opts, 'H', nan(n*m, 1));
% Allocate space for required number of handles
if length(H)<n*m, H = [H; nan(n*m-length(H),1)]; end

% Error if plot number is outside range of subplots
if p>n*m, error('Index exceeds number of subplots'); end;

% Update position of all existing subplots
for PlotNum = 1:n*m
    if ishandle(H(PlotNum))
        R = computePos(PlotNum);
        set(H(PlotNum), 'Position', R);
    end
end

% If the current subplot does not exist yet
if ~ishandle(H(p))
    % Compute its position and create it
    R = computePos(p);
    H(p) = axes('position', R);

else
% If the current subplot already exists 
    % Activate it
    ch = get(gcf, 'Children');
    ch(ch==H(p)) = [];
    set(gcf, 'Children', [H(p);ch]);
    set(gcf, 'CurrentAxes', H(p));
end

opts.H = H;
% Save parameters as userdata
set(gcf, 'Userdata', opts);

%--------------------------------------------------------------------------
function R = computePos(p)

    % Grid indices for the subplot
    [pn,pm] = ind2sub([n,m], p);

    % Compute height and width of subplot
    am = (1-BorderTop-BorderBottom-(m-1)*SpacingY)*Heights/sum(Heights);
    an = (1-BorderLeft-BorderRight-(n-1)*SpacingX)*Widths/sum(Widths);

    % Compute coordinates of subplot
    R = zeros(1,4);
    R(1) = BorderLeft+(pn-1)*(SpacingX)+sum(an(1:pn-1));
    R(2) = BorderBottom+(m-pm)*(SpacingY)+sum(am(pm+1:end));
    R(3) = an(pn);
    R(4) = am(pm);
end

end


function out = mgetopt(varargin)
% MGETOPT Parser for optional arguments
% 
% Usage
%   Get a parameter structure from 'varargin'
%     opts = mgetopt(varargin);
%
%   Get and parse a parameter:
%     var = mgetopt(opts, varname, default);
%        opts:    parameter structure
%        varname: name of variable
%        default: default value if variable is not set
%
%     var = mgetopt(opts, varname, default, command, argument);
%        command, argument:
%          String in set:
%          'instrset', {'str1', 'str2', ... }
%
% Example
%    function y = myfun(x, varargin)
%    ...
%    opts = mgetopt(varargin);
%    parm1 = mgetopt(opts, 'parm1', 0)
%    ...

% Copyright 2007 Mikkel N. Schmidt, ms@it.dk, www.mikkelschmidt.dk

if nargin==1
    if isempty(varargin{1})
        out = struct;
    elseif isstruct(varargin{1})
        out = varargin{1}{:};
    elseif isstruct(varargin{1}{1})
        out = varargin{1}{1};
    else
        out = cell2struct(varargin{1}(2:2:end),varargin{1}(1:2:end),2);
    end
elseif nargin>=3
    opts = varargin{1};
    varname = varargin{2};
    default = varargin{3};
    validation = varargin(4:end);
    if isfield(opts, varname)
        out = opts.(varname);
    else
        out = default;
    end
    
    for narg = 1:2:length(validation)
        cmd = validation{narg};
        arg = validation{narg+1};
        switch cmd
            case 'instrset',
                if ~any(strcmp(arg, out))
                    fprintf(['Wrong argument %s = ''%s'' - ', ...
                        'Using default : %s = ''%s''\n'], ...
                        varname, out, varname, default);
                    out = default;
                end
            case 'dim'
                if ~all(size(out)==arg)
                    fprintf(['Wrong argument dimension: %s - ', ...
                        'Using default.\n'], ...
                        varname);
                    out = default;
                end
            otherwise,
                error('Wrong option: %s.', cmd);
        end
    end
end
end