function meqplotQ(L,varargin)
%MEQPLOTQ: compare profiles and geometry for static equilibria with the same L structure
% meqplotQ(L,LY1,LY2)
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

% post processing was run for at least one LY
ipost = any(cellfun(@(x)isfield(x,'iqQ'),varargin)); 

if ipost, nr = 4; else, nr=2; end; % rows for plots
hs(1) = subplot(1 ,3, 1);
hs(2) = subplot(nr,3, 2);
hs(3) = subplot(nr,3, 3);
hs(4) = subplot(nr,3, 5);
hs(5) = subplot(nr,3, 6);
if ipost
  hs(6) = subplot(nr,3, 8);
  hs(7) = subplot(nr,3, 9);
  hs(8) = subplot(nr,3,11);
  hs(9) = subplot(nr,3,12);
end

meqgplot(L.G,hs(1),'vl'); hold(hs(1),'on');

for ii=1:numel(varargin)
  LY = varargin{ii};
  
  % postprocessing is available for this LY
  ipost = isfield(LY,'iqQ'); 
  
  %% flux plots
  F = linspace(LY.FA,2*LY.FB-LY.FA,21);
  
  ax=hs(1);
  col = ax.ColorOrder(ii,:);
  
  contour(ax,L.G.rx,L.G.zx,LY.Fx,F,'color',col);
  contour(ax,L.G.rx,L.G.zx,LY.Fx,[1 1]*LY.FB,'color',col,'linewidth',2);
  
  plot(ax,LY.rA,LY.zA,'ow');
  plot(ax,LY.rB,LY.zB,'gs');
  
  axis(ax,'equal','tight');
  
  %% profiles
  ax = hs(2);
  plot(ax,L.pQ.^2,mu0*LY.PpQ,'color',col);
  hold(ax,'on');
  ylabel(ax,'\mu_0p'''); xlabel(ax,'\psi_N');
  
  ax = hs(3);
  plot(ax,L.pQ.^2,LY.PQ/1e3,'color',col);
  hold(ax,'on');
  ylabel(ax,'p [kPa]'); xlabel(ax,'\psi_N');
  
  ax = hs(4);
  plot(ax,L.pQ.^2,LY.TTpQ,'color',col);
  ylabel(ax,'TT''');xlabel(ax,'\psi_N');
  hold(ax,'on');
  
  ax = hs(5);
  plot(ax,L.pQ.^2,LY.TQ,'color',col);
  ylabel(ax,'T [Wb]');xlabel(ax,'\psi_N');
  hold(ax,'on');
  
  if ipost
    ax = hs(6);
    plot(ax,L.pQ.^2,1./LY.iqQ,'color',col);
    ylabel(ax,'q'); xlabel(ax,'\psi_N');
    hold(ax,'on');
    
    
    ax = hs(7);
    plot(ax,L.pQ.^2,LY.kappa,'color',col);
    ylabel(ax,'\kappa'); xlabel(ax,'\psi_N');
    hold(ax,'on');
    
    ax = hs(8);
    plot(ax,L.pQ.^2,LY.jtorQ*1e-6,'color',col);
    ylabel(ax,'j_{tor} [MA/m^2]');  xlabel(ax,'\psi_N');
    hold(ax,'on');
    
    ax = hs(9);
    plot(ax,L.pQ.^2,LY.delta,'color',col);
    ylabel(ax,'\delta');  xlabel(ax,'\psi_N');
    hold(ax,'on');
  end
end
