function hax = meqplotfull(L,LY,varargin)
% Plot time-evolution with R,Z,Ip
% and time-evolving coil currents, voltages and equations, with limits
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

p=inputParser;
p.addParameter('hax',[]);
parse(p,varargin{:}); P = p.Results;

if isempty(P.hax)
  hax(1) = subplot(1,2,1);
  hax(2) = subplot(6,2,2);
  hax(3) = subplot(6,2,4);
  hax(4) = subplot(6,2,6);
  hax(5) = subplot(6,2,8);
  hax(6) = subplot(6,2,10);
  hax(7) = subplot(6,2,12);
  
  thax = hax(5:7); % axes with time
  redraw = true;
else
  % Find axes handles to replace data
  hax = P.hax;
  hVbar = findobj(hax(2),'Tag','Vbar');
  hIbar = findobj(hax(3),'Tag','Ibar');
  hPbar = findobj(hax(4),'Tag','Pbar');
  hTbar = findobj(hax(4),'Tag','Tbar');
  hzA = flipud(findobj(hax(5),'Tag','zA'));
  hrA = flipud(findobj(hax(6),'Tag','rA'));
  hIp = flipud(findobj(hax(7),'Tag','Ip'));
  redraw = false;
end

%% meqplotfancy
meqplotfancy(L,LY,'parent',hax(1));

%% Data
if isfield(L.G,'Ihig'), np = numel(L.G.Ihig); else, np = 0; end
iTa = ~isinf(L.G.Talim); % Coils for which thermal limits exist
nt = sum(iTa);

% defaults
if isfield(LY,'Va'), Varel = LY.Va./L.G.Valim; else; Varel = zeros(size(LY.Ia)); end
if ~isfield(LY,'Iarel'), LY.Iarel = LY.Ia./L.G.Ialim; end
if ~isfield(LY,'Parel'), LY.Parel = zeros(np,numel(LY.t)); end
if ~isfield(LY,'Tarel'), LY.Tarel = zeros(numel(L.G.Talim),numel(LY.t));end

%%  plots
if redraw
  [tick,lims,label] = coilticks(L.G.dima);
  %% Va plot
  ax = hax(2);
  hold(ax,'on');
  coilpatch(ax,lims);
  bar(1:L.G.na,Varel,'parent',ax,'FaceColor','b',...
    'Tag','Vbar');
  
  stairs(ax,-0.5+(1:L.G.na+1), ones(1,L.G.na+1),'k')
  stairs(ax,-0.5+(1:L.G.na+1),-ones(1,L.G.na+1),'k')
  
  ylabel(ax,'V_a/V_{a,lim}');
  set(ax,'Ylim',1.1*[-1 1],'Xlim',0.5+[1,L.G.na])
  set(ax,'Xtick',tick,'XtickLabel',label)
  box(ax,'on');
  
  
  %% Ia plot
  ax = hax(3);
  hold(ax,'on');
  coilpatch(ax,lims);
  bar(1:L.G.na,LY.Iarel,'parent',ax,'FaceColor','b',...
    'Tag','Ibar');
  
  stairs(ax,-0.5+(1:L.G.na+1), ones(1,L.G.na+1),'k')
  stairs(ax,-0.5+(1:L.G.na+1),-ones(1,L.G.na+1),'k')
  
  ylabel(ax,'I_a/I_{a,lim}');
  set(ax,'Ylim',1.1*[-1 1],'Xlim',0.5+[1,L.G.na])
  set(ax,'Xtick',tick,'XtickLabel',label)
  box(ax,'on');
  
  
  %% Coil protection limit plots
  ax = hax(4);
  hold(ax,'on');
  
  if np>0
    coilpatch(ax,[0 np np+nt]);
    
    bar(1:np,LY.Parel,'parent',ax,'FaceColor','b',...
      'Tag','Pbar');
    
    Iprotmin = min(L.G.Ihig,-L.G.Ilow);
    hi = L.G.Ihig./Iprotmin;
    lo = L.G.Ilow./Iprotmin;
    stairs(ax,-0.5+(1:np+1),[hi;hi(end)],'k')
    stairs(ax,-0.5+(1:np+1),[lo;lo(end)],'k')
  end
  % Thermal protection limits
  
  if nt>0
    bar(np+(1:nt),LY.Tarel(~isinf(L.G.Talim)),...
      'parent',ax,'Tag','Tbar');
    
    stairs(ax,-0.5+np+(1:nt+1), [1 1],'k')
    stairs(ax,-0.5+np+(1:nt+1),-[1 1],'k')
  end
  
  if np+nt > 0
    ylabel(ax,'Protections');
    set(ax,'Ylim',[-1.1 1.1],'Xlim',-0.5+[1,np+nt+1])
    box(ax,'on');
    iplabel = 1:2:np;
    protlabel = [cellstr(num2str((iplabel)','%d'));'';L.G.dima{iTa}(1)];
    set(ax,'Xtick',[iplabel,np+(1:nt)],'XtickLabel',protlabel);
  end
  %% R,Z,Ip plots
  ax = hax(5);
  plot(ax,LY.t,LY.zA,'.-','Tag','zA');
  hold(ax,'on');
  ylabel(ax,'z_A')
  
  ax = hax(6);
  plot(ax,LY.t,LY.rA,'.-','Tag','rA');
  hold(ax,'on');
  ylabel(ax,'r_A')
  
  ax = hax(7);
  plot(ax,LY.t,LY.IpD/1e3,'.-','Tag','Ip');
  hold(ax,'on');  axis(ax,'tight')
  xlabel(ax,'t[s]'); ylabel(ax,'I_{pD}[kA]')
  linkaxes(thax,'x')
  set(thax,'Xlim',LY.t+[0 1e-3]);
  
else
  % update bar data
  set(hVbar,'YData',Varel);
  set(hIbar,'YData',LY.Iarel);
  set(hPbar,'YData',LY.Parel);
  set(hTbar,'YData',LY.Tarel);
  
  % update time data
  for iA=1:LY.nA
    t = [get(hzA(iA),'Xdata'),LY.t]; % extend time
    set(hzA(iA), 'Xdata',t, 'YData',[get(hzA(iA),'Ydata'),LY.zA(iA)]);
    set(hrA(iA), 'Xdata',t, 'YData',[get(hrA(iA),'Ydata'),LY.rA(iA)]);
    set(hIp(iA), 'Xdata',t, 'YData',[get(hIp(iA),'Ydata'),LY.IpD(iA)/1e3]);
  end
  for iB = LY.nA+1:L.nD
    set(hIp(iB), 'Xdata',t, 'YData',[get(hIp(iB),'Ydata'),LY.IpD(iB)/1e3]);
  end
  
  hga = get(hzA(1),'parent'); % handle of axis containing t plot
  tlim =   get(hga,'Xlim');
  if LY.t > tlim(2)
    tlim(2) = LY.t+0.5*(tlim(2)-tlim(1));
    set(hga,'Xlim',tlim); % extend time range;
  end
end
end

function [tick,lims,coiltypes] = coilticks(dima)
coiltypes = unique(extractBefore(string(dima),'_'));
tick = zeros(numel(coiltypes),1);
lims = zeros(numel(coiltypes)+1,1);
lims(1) = 0; % init
for icoiltype = 1:numel(coiltypes)
  icoil = find(startsWith(dima,coiltypes{icoiltype}));
  lims(icoiltype+1) = icoil(end);
  tick(icoiltype) = mean(icoil);
end
end

function coilpatch(ax,lims)
gy = 0.8*[1 1 1]; lgy = 0.9*[1 1 1];
for ii=1:numel(lims)-1
  xx = 0.5 + [lims(ii), lims(ii+1) lims(ii+1) lims(ii)];
  yy = 0.99 * [-1 -1 1 1];
  if isodd(ii),col=gy; else; col=lgy; end
  hp = patch(xx,yy,col,'parent',ax,'EdgeColor','none');
end
end