function meqcompare(L,varargin)
%MEQCOMPARE compare static equilibria with the same L structure
% example:
%    meqcompare(L,LY1,LY2)
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

 hs(1) = subplot(1,2,1);
 hs(2) = subplot(6,2,2);
 hs(3) = subplot(6,2,4);
 hs(4) = subplot(6,2,6);
 hs(5) = subplot(6,2,8);
 hs(6) = subplot(6,2,10);
 hs(7) = subplot(6,2,12);

 meqgplot(L.G,hs(1),'vlw'); hold(hs(1),'on');
 
 nLY = numel(varargin);
 wbar=0.8/nLY; % bar graph widths

 assert(nLY<5,'compare a maximum of 4 LY');
 LineStyleOrder = {'-','--','-.',':'};
 for ii=1:nLY
  LY = varargin{ii};
  assert(numel(LY.t)==1,'Can only compare scalar LY structures, use meqxk() to slice')
  %% flux plots
  ax=hs(1); hold(ax,'on');
  col = ax.ColorOrder(ii,:);
  sty = LineStyleOrder{ii};
  
  FFB = repmat(LY.FB,1,max(3-numel(LY.FB),1));
  if ii==1; F=linspace(min(LY.Fx(:)),max(LY.Fx(:)),21); end; % same contour values for all cases
    
  contour(ax,L.G.rx,L.G.zx,LY.Fx,F,'color',col,'LineStyle',sty);
  contour(ax,L.G.rx,L.G.zx,LY.Fx,FFB,'color',col,'linewidth',2,'LineStyle',sty);
  
  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,mu0*LY.PpQ,'color',col,'LineStyle',sty);
  hold(ax,'on');
  title(ax,'\mu_0p''')
  
  ax = hs(3);
  plot(ax,L.pQ,LY.TTpQ,'color',col,'LineStyle',sty);
  title(ax,'TT''')
  hold(ax,'on');
  
  %% Currents
  ax = hs(4);
  bar(ax,(1:L.G.na)+(ii-1)*wbar+0.1,LY.Ia,wbar,'FaceColor',col)
  hold(ax,'on');
  title(ax,'Ia');

  ax = hs(5);
  plot(ax,1:L.G.nv,LY.Iv,'Color',col,'LineStyle',sty)
  hold(ax,'on');
  title(ax,'Iv');

  
  %% measurements
  ax = hs(6);
  bar(ax,(1:L.G.nm)+(ii-1)*wbar+0.1,LY.Bm,wbar,'FaceColor',col)
  hold(ax,'on');
  title(ax,'Bm');

  ax = hs(7);
  bar(ax,(1:L.G.nf)+(ii-1)*wbar+0.1,LY.Ff,wbar,'FaceColor',col)
  hold(ax,'on');
  title(ax,'Ff');
 end