function hax=meqplott(L,LY)
%MEQPLOTT - A simple plot of a single time slice
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

% plot a single timeslice
assert(~isempty(LY),'cannot plot empty LY');
if isfield(LY,'t'), t = LY.t; else t=NaN; end
if isfield(LY,'niter'), niter = LY.niter; else niter = NaN; end
if isfield(LY,'shot')
  shot = LY.shot;
elseif isfield(L.P,'shot')
  shot = L.P.shot;
else
  shot = NaN;
end

assert(numel(LY.t)==1,'can plot only for single timeslice, use meqxk() to slice');
Iy = LY.Iy;
Fx = LY.Fx;
rX = LY.rX; zX = LY.zX;

if isfield(LY,'Iy') && any(LY.Iy(:))
  Iy = LY.Iy;
  FA = LY.FA; FB = LY.FB;
  rA = LY.rA; zA = LY.zA;
  rB = LY.rB; zB = LY.zB;
  rX = LY.rX; zX = LY.zX;
else
  [FA,FB,rA,rB,zA,zB,rX,zX] = deal(nan);
end

if isfield(LY,'dz'); dz=LY.dz; else, dz=zeros(size(rA)); end
if isfield(L,'bh') 
  rbh=L.bh([1 3 3 1 1]);
  zbh=L.bh([2 2 4 4 2]);
else
  rbh=[]; zbh=[];
end
if isfield(L,'rh'); rh=L.rh; else, rh=[]; end
if isfield(L,'zh'); zh=L.zh; else, zh=[]; end

rl = [L.G.rl;L.G.rl(1)]; zl = [L.G.zl;L.G.zl(1)];
% Symbols-lines to plot on the cross section
mvis = {'w-', 'r.', 'c-o', 'om','om','.-c',  '-y',   'oy','wx'};
rvis = {[]     []    []     rB   rA   rl     rbh      rh   rX };
zvis = {[]     []    []     zB   zA   zl     zbh      zh   zX };

if ~isempty(LY.aq)
  aq = LY.aq; rq=[];zq=[];
  %% contours
  for iD = 1:LY.nA
    rqD = rA(iD) + aq(:,:,iD).*L.crq;
    zqD = zA(iD) + aq(:,:,iD).*L.czq;
    rq = [rq,rqD([1:end 1],:)];
    zq = [zq,zqD([1:end 1],:)];
  end
  rvis{2} = rq;
  zvis{2} = zq;
end

if ~isempty(LY.aW)
  aW = LY.aW;
  rvis{3} = L.G.rW' + [zeros(size(aW)) -aW.*cos(L.G.oW)]';
  zvis{3} = L.G.zW' + [zeros(size(aW)) +aW.*sin(L.G.oW)]';
end

if ~isfield(LY,'Ip'); LY.Ip = sum(Iy(:)); end
if ~isfield(LY,'bp'); LY.bp = NaN; end
if ~isfield(LY,'li'); LY.li = NaN; end

if niter == 1
  % on first iteration plot the contour lines of Iy
  [rvis{1},zvis{1}] = mycontourc(L.ry,L.zy,double(abs(Iy)),linspace(.1,.95,5)*max(double(abs(Iy(:)))));
else
  rvis{1} = [];
  zvis{1} = [];
end

% Optional plots
if ~isfinite(FB) % if meqpdom could not find a valid LCFS point
  F = 21;
elseif ~isempty(FA) % plot specific contour lines depending on FB,FBA
  sIp = sign(LY.Ip);
  FA_ = sIp*max(sIp*FA);
  FB_ = sIp*min(sIp*FB);
  FBA = (FB_ - FA_);
  F = FB_-[linspace(-1,0,11), L.fq(end:-1:1),linspace(1,2.5,16) ]*FBA;
else 
  F = 31;
end
  
rztext = cell(1,LY.nA);
dz(end+1:LY.nA) = 0;
for iA = 1:LY.nA
  rztext{iA} = sprintf('r,z=%5.3f,%+6.3fm%+4.1fmm:%d',rA(iA),zA(iA),dz(iA)*1000,sum(LY.Opy(:) == iA));
end

meqvis([{sprintf('%s#%d %6.4fs/%d',L.P.tokamak,shot,t,niter)} ...
  rztext ...
  {sprintf('Ip=%.0fkA bp=%.2f li=%.2f',LY.Ip*1e-3,LY.bp,LY.li)}], ...
  L.G.rx,L.G.zx,Fx,F,FB,...
  rvis,zvis,mvis,...
  {})

if nargout>0
  hax=gca;
end

end
