function varargout = meqplotevo(L,varargin)
%MEQPLOTEVO Plot time evolution of meq LY structures
% meqplotevo(L,LY1,LY2,LY3,'parameter',value)
% Optional parameters inputs:
%   argplot: Cell array with LY fields to plot
%   legstr:  Cell array with strings for legend
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

varargins = {};

assert(nargin>1,'insufficient arguments')
for ii=1:numel(varargin)
  if isstruct(varargin{ii})
    LYs{ii} = varargin{ii}; %#ok<AGROW>
  else
    varargins = varargin(ii:end);
    break
  end
end

if L.P.idoublet, Ip='IpD'; else, Ip='Ip'; end

iP = inputParser; % initialize parser
addParameter(iP,'argplot',{'Ia','Va',Ip,'zA','rA','Wk','bp','bt','li','qA'})
addParameter(iP,'legstr',{''});
parse(iP,varargins{:}); % parse the inputs, input xyz is now available in P.xyz
P = iP.Results;

% plot and compare evolution of time-varying meq LY results
np = numel(P.argplot); % number of plots
nc = 2; nr = ceil(np/nc); % rows/columns of plots

clf; hax = zeros(nr,nc);
for ii = 1 : np
  hax(ii) = subplot(nr,nc,ii);
end
nLY = numel(LYs);


cols = get(gca,'ColorOrder');

for iLY = 1:nLY
  LY = LYs{iLY};
  col = cols(iLY,:);
 
  for iplot=1:np
    ax = hax(iplot);
    
    myplot = P.argplot{iplot};
 
    if isfield(LY,myplot)
      t = LY.t;
      data = LY.(myplot);
      if size(data,2)==1, t=1; end
      plot(ax,t,data,'Color',col);
      hold(ax,'on');
    end
    
    title(ax,myplot)
  end
end
linkaxes(hax,'x');
set(hax(1),'XLim',[t(1),t(end)]);
if ~isempty(P.legstr)
  legend(P.legstr)
end

if nargout==1
  varargout=hax;
end
end
