function meqplotfield(L,LY,varargin)
%MEQPLOTFIELD Plot magnetic field from LY structure
% function meqplotfield(L,LY,varargin)
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

p=inputParser; p.KeepUnmatched = true; 

p.addParameter('parent',gca,@(x) ishandle(x));
p.addParameter('BpLineStyle','--'); % line style for Bp
p.addParameter('BpColor','r'); % color for Bp
p.addParameter('BpLevels',[0:1e-3:5e-3]); % levels to plot
p.addParameter('BpLabels',true); % plot contour labels
p.addParameter('BpScale',1e3); % scale
p.addParameter('plotBp',true); % add |Bp| contour plot 
p.addParameter('plotArrows',true);
p.addParameter('vacuum',false); % plots vacuum field Br0x,Bz0x is available
p.addParameter('izgrid',false); % plot on extended grid if present
parse(p,varargin{:}); P = p.Results;

hax = P.parent;

assert(numel(LY.t)==1,'can only plot for one timeslice, use meqxk() to slice')

if P.izgrid
  assert(all(isfield(LY,{'Brz','Bzz','Fz'})),'must run with izgrid=true to plot fields on extended grid');
  rx=L.rz; zx=L.zz; 
  if P.vacuum && LY.Ip~=0
    Brx = LY.Br0z; Bzx = LY.Bz0z; Fx  = LY.F0z;
  else
    Brx = LY.Brz; Bzx = LY.Bzz; Fx  = LY.Fz;
  end   
else % default x grid
  rx=L.rx; zx=L.zx; 
  if P.vacuum && LY.Ip~=0
    assert(all(isfield(LY,{'Br0x','Bz0x','F0x'})),'missing fields Br0x,Bz0x or F0x, required to plot vacuum field if Ip~=0')
    Brx = LY.Br0x; Bzx = LY.Bz0x; Fx  = LY.F0x; 
  else
    Brx = LY.Brx; Bzx = LY.Bzx; Fx  = LY.Fx;
  end
end  

axis(hax,'equal','tight');
hold(hax,'on');
if P.plotBp
  [c,h] = contour(rx,zx,sqrt(Brx.^2+Bzx.^2)*P.BpScale,P.BpLevels*P.BpScale,...
    'LineStyle',P.BpLineStyle,'Color',P.BpColor,'parent',hax);
  if P.BpLabels
    clabel(c,h);
  end
end

if P.plotArrows
  meqplotarrows(rx,zx,Fx,Brx,Bzx,varargin{:},'plotContour',true,'parent',hax);
end

end
