function meqplotD(L,LY,varargin)
% meqplotD(L,LY,LY2..)
% Plot droplet evolution
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

LYs{1} = LY; % in cell
for ii=1:numel(varargin)
  if isstruct(varargin{ii})
    LYs{1+ii} = varargin{ii};
  else
    %varargin = varargin(ii:nargin); break; % later opts
  end
end

cols = get(gca,'ColorOrder');
for iLY = 1:numel(LYs)
  col = cols(iLY,:);

  LY = LYs{iLY};
  % Totals
  subplot(221); plot(LY.t,LY.Ip/1e3); title('Ip [kA]'); hold on;
  
  % Per domain
  for iD=1:size(LY.IpD,1)
    switch iD
      case 1, sty='-' ;
      case 2, sty='--';
    end
    
    m = LY.lB(iD,:);
    
    subplot(222); plot(LY.t(m),LY.IpD(iD,m)/1e3,'LineStyle',sty,'color',col); hold on;
    title('IpD [kA]')
    
    subplot(223); plot(LY.t(m),LY.rA(iD,m),'LineStyle',sty,'color',col); hold on;
    title('r'); ylim([min(L.G.rv),max(L.G.rv)]);
    
    subplot(224); plot(LY.t(m),LY.zA(iD,m),'LineStyle',sty,'color',col); hold on;
    title('z'); ylim([min(L.G.zv),max(L.G.zv)]);
    
  end
end
end