function varargout=fbtxplot(L,LX,varargin)
% function hax = fbtxplot(L,LX,varargin)
%
% Plot geometry constraints in LX
%
% Optional parameter-value pairs
%   'parent'     : axis handle on which to plot (default:gca)
%   'color'      : color for this plot
%   'legend'     : show legend or not (default:false)
%   'legendopts' : options for legend
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

p=inputParser;
p.addParameter('parent',gca,@(x) ishandle(x));
p.addParameter('color','b',@(x) ischar(x)||isnumeric(x));
p.addParameter('legend',true);
p.addParameter('legendopts',{});

parse(p,varargin{:}); P = p.Results;

hax=P.parent;

assert(size(LX.gpr,2)==1,'fbtxplot accepts only one timeslice, use meqxk() to slice')
hp=gobjects(0);

%%
meqgplot(L.G,hax,'lv'); hold on;

for it=1:numel(LX.t)
  hp = fbtxplot_oneslice(hax,hp,L,meqxk(LX,it),P);
end

if P.legend
  legend(hp,'location','southoutside','autoupdate','off',P.legendopts{:});
end

if nargout>0
  varargout{1} = hax;
end
end


function hp = fbtxplot_oneslice(hax,hp,L,LX,P,varargin)
nS = size(LX.gpr,1);

iWf = LX.gpfd.*LX.gpfe;
iWb = LX.gpbd.*LX.gpbe;
iWc = LX.gpcd.*LX.gpce;
iWv = LX.gpvd.*LX.gpve;

for iS = 1:nS
  % Loop over constraints and add marker properties
  
  MarkerFull = [];
  Marker = '';
  MarkerSize = [];
  DisplayName = {};
  
  rr = LX.gpr(iS); zz=LX.gpz(iS);
  
  if ~any(isnan([iWf(iS),LX.gpfa(iS),LX.gpfb(iS)]))  % flux constraint
    MarkerFull(end+1) = LX.gpb(iS); %#ok<*AGROW>
    if LX.gpfe(iS)~=0
      if LX.gpfb(iS)~=0, Marker(end+1) = 'o'; ms=5; DisplayName{end+1} = 'Relative flux cost';
      else,              Marker(end+1) = '<'; ms=5; DisplayName{end+1} = 'Absolute flux cost';
      end
    else
      if LX.gpfb(iS)~=0, Marker(end+1) = 's'; ms=8; DisplayName{end+1} = 'Relative flux constraint';
      else,              Marker(end+1) = '>'; ms=8; DisplayName{end+1} = 'Absolute flux constraint';
      end
    end
    MarkerSize(end+1) = ms;
  end
  
  isBr = ~any(isnan([iWb(iS),LX.gpbr(iS)])); % Br constraint
  isBz = ~any(isnan([iWb(iS),LX.gpbz(iS)])); % Bz constraint
  isBa = ~any(isnan([iWb(iS),LX.gpba(iS)])); % Bz constraint
  
  if isBr && isBz && LX.gpbr(iS)==0 && LX.gpbz(iS)==0
    % X point
    Marker(end+1)     = 'x';
    MarkerSize(end+1) = 15 ;
    MarkerFull(end+1) = false;
    DisplayName{end+1} = 'X point';
  elseif isBr || isBz  || isBa % one of Br or Bz or Ba constrained
    Marker(end+1) = '+';
    DisplayName{end+1} = sprintf('B field %s',ifc(LX.gpbe==0,'constraint','cost'));
    MarkerFull(end+1) = false;
    MarkerSize(end+1) = 15 ;
  end
  
  if ~isnan(iWc(iS)) && any(~isnan([LX.gpcr(iS),LX.gpcz(iS),LX.gpca(iS)])) % Hessian constraint
    Marker(end+1) = '*';
    DisplayName{end+1} = sprintf('Flux Hessian %s',ifc(LX.gpce==0,'constraint','cost'));
    MarkerFull(end+1)  = false;
    MarkerSize(end+1)  = 15 ;
  end
  
  if ~isnan(iWv(iS)) && any(~isnan([LX.gpvrr(iS),LX.gpvrz(iS),LX.gpvzz(iS)]))  % VAcuum field constraint
    Marker(end+1) = 'd';
    DisplayName{end+1} = sprintf('Vacuum field %s',ifc(LX.gpve==0,'constraint','cost'));
    MarkerFull        = false;
    MarkerSize(end+1) = 10 ;
  end
  
  hp = doplot(hax,hp,rr,zz,P.color,Marker,MarkerFull,MarkerSize,DisplayName,varargin{:});
end
end

function hp = doplot(hax,hp,rr,zz,Color,Marker,MarkerFull,MarkerSize,DisplayName,varargin)
% plot all the accumulated markers with their properties
for ii=1:numel(Marker)
  if MarkerFull(ii), MarkerFaceColor = Color; else, MarkerFaceColor='none'; end
  hn=plot(hax,rr,zz,'Marker',Marker(ii),'MarkerEdgeColor',Color,'Markersize',MarkerSize(ii),...
    'MarkerFaceColor',MarkerFaceColor,'DisplayName',DisplayName{ii},varargin{:});
  if isempty(hp) || ~contains(hn.DisplayName,{hp.DisplayName})
    hp(end+1)=hn;
  end
end
end
