function varargout = fbtxdisp(L,LX,selection)
% function s = fbtxdisp(L,LX,selection)
%
% Textual display of FBT constraints.
% If called without output argument, displays in the command window.
% If passing output argument s, returns cell array with strings
%
% selection: optional argument to choose what to display among the following categories:
%            'current','dipole','limits','geometry'. By default all are selected.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

assert(numel(LX.t)==1,'supported for 1 timeslice only, use meqxk() to slice')
if nargout>0, s={}; else, s=[]; end % init cell array for storing messages
ind0 = 25; % indentation
ind1 = ind0-18; % smaller indentation

selection_all = {'currents','dipole','limits','geometry'};
if nargin<3, selection = selection_all; % default
else
  if ~iscell(selection), selection={selection}; end
  assert(all(cellfun(@(x) contains(x,selection_all),selection)),'invalid selection, allowed: %s',selection_all)
end
if isempty(selection); return; end


%% Header
s=vprintf(s,'\nFBTX cost function / constraints display\n');
s=vprintf(s,'  Tokamak: %s, shot#%d, t=%3.3f\n',L.P.tokamak,LX.shot(1),LX.t(1));

%% Inverse weights
iWa = LX.gpid.*LX.gpie;
iWd = LX.gpdd.*LX.gpde;
iWf = LX.gpfd.*LX.gpfe;
iWB = LX.gpbd.*LX.gpbe;
iWC = LX.gpcd.*LX.gpce;
iWV = LX.gpvd.*LX.gpve;

%% current constraints

if contains('currents',selection)
  s=vprintf(s,'\nCoil current cost / equality constraints:\n');
  
  isela = find(isfinite(iWa));
  % cost function or equality constraint per coil current
  for ia = isela'
    s = vprintf(s,'%-18s Current:  Ia = %2.2e[A]                               1/w=%2.2e\n',L.G.dima{ia},LX.gpia(ia),iWa(ia));
  end
  if isempty(ia); s=vprintf(s,'   none\n'); end
end

if contains('dipoles',selection)
  s=vprintf(s,'\nDipole cost / equality constraints:\n');
  % Dipoles
  for id = find(isfinite(iWd))'
    s = vprintf(s,'%-18s Dipole :  (%4.2g*Ia(%2d)-%4.2g*Ia(%2d))=0      1/w=%2.2e\n',...
      sprintf('%s-%s',L.G.dima{id+1},L.G.dima{id}),LX.gpdw(id+1),id+1,LX.gpdw(id),id,LX.gpdd*LX.gpde(id));
  end
  if isempty(id); s=vprintf(s,'   none\n'); end
end

if contains('limits',selection)
  s=vprintf(s,'\nCoil current limits:\n');
  if all(isinf([L.P.limu(:);L.P.liml(:)]))
    s=vprintf(s,'   none\n');
  else
    % Inequality constraints: liml*limm < limc*Ia < limu*limm
    limu = L.P.limu.*L.P.limm; % apply margin
    liml = L.P.liml.*L.P.limm;
    for ic = 1:size(L.P.limc,1)
      limc = L.P.limc(ic,:);
      if all(isinf([limu(ic),liml(ic)]))||~any(limc), continue; end
      
      % for each row, find the coils involved
      iinvolved = find(limc~=0);
      expr = '';
      for iinv = iinvolved
        expr = sprintf('%s %+5.2g*%s',expr,limc(iinv),L.G.dima{iinv});
      end
      if ~isinf(liml(ic))
        expr = sprintf('%4.2e <= %s',liml(ic),expr);    % lower limit
      end
      if ~isinf(limu(ic))
        expr = sprintf('%-30s <= %4.2g',expr,limu(ic)); % upper limit
      end
      s = vprintf(s,'%s\n',expr); % print
    end
  end
end

%% display per shape point
if contains('geometry',selection)
  
  rr=NaN; zz=NaN; % init
  
  s=vprintf(s,'\nMagnetic geometry cost / equality constraints:\n');
  nS = numel(LX.gpr); %number of points
  for iS=1:nS
    if any(isnan([LX.gpr(iS),LX.gpz(iS)])) % not a valid R,Z point
      continue;
    end
    
    if ~(LX.gpr(iS)==rr && LX.gpz(iS)==zz) % point not seen yet
      rr = LX.gpr(iS); zz=LX.gpz(iS); % track previously displayed point
      s = vprintf(s,'r=%+4.3f z=%+4.3f ',LX.gpr(iS),LX.gpz(iS)); % display R,Z
      nx = false; % reset
    end
    
    if ~any(isnan([iWf(iS),LX.gpfa(iS),LX.gpfb(iS)])) % flux constraint
      if ~nx, ind = ind1; nx=true; else, ind=ind0; end
      if LX.gpb(iS), sB=sprintf(' B ');  else,sB=sprintf('   '); end    % indicate that point is on LCFS
      s=vprintf(s,'%*s:   MSa*Ia +  MSy*Iy = %5.2g + %4.2g*Fb                 1/w=%2.2e\n',ind,[sB,'Psi'],LX.gpfa(iS),LX.gpfb(iS),iWf(iS));
    end
    
    % BR,BZ constraints
    if ~any(isnan([iWB(iS),LX.gpbr(iS)]))
      if ~nx, ind = ind1; nx=true; else, ind=ind0; end
      s=vprintf(s,'%*s:  BrSa*Ia + BrSy*Iy = %5.2g                           1/w=%2.2e\n',ind,'Br',LX.gpbr(iS),            iWB(iS));
    end
    if ~any(isnan([iWB(iS),LX.gpbz(iS)]))
      if ~nx, ind = ind1; nx=true; else, ind=ind0; end
      s=vprintf(s,'%*s:  BzSa*Ia + BzSy*Iy = %5.2g                           1/w=%2.2e\n',ind,'Bz',LX.gpbz(iS),            iWB(iS));
    end
    
    % strike point angle constraint
    if ~any(isnan([iWB(iS),LX.gpba(iS)]))
      if ~nx, ind = ind1; nx=true; else, ind=ind0; end
      s=vprintf(s,'%*s: %4.1f*(BrSa*Ia+BrSy*Iy) + %4.1f*(BzSa*Ia+BzSy*Iy) = 0  1/w=%2.2e\n',ind,'BN',sin(LX.gpba(iS)),cos(LX.gpba(iS)),iWB(iS));
    end
    
    % Second order flux derivative constraint
    if ~any(isnan([iWC(iS),LX.gpcr(iS)]))
      if ~nx, ind = ind1; nx=true; else, ind=ind0; end
      s=vprintf(s,'%*s: MrrSa*Ia +MrrSy*Iy = %5.2g                           1/w=%2.2e\n',ind,'d2r psi',LX.gpcr(iS),iWC(iS));
    end
    if ~any(isnan([iWC(iS),LX.gpcz(iS)]))
      if ~nx, ind = ind1; nx=true; else, ind=ind0; end
      s=vprintf(s,'%*s: MrzSa*Ia +MrzSy*Iy = %5.2g                           1/w=%2.2e\n',ind,'drdz psi',LX.gpcz(iS),iWC(iS));
    end
    
    % Hessian angle constraint
    if ~any(isnan([iWC(iS),LX.gpca(iS)]))
      if ~nx, ind = ind1; nx=true; else, ind=ind0; end
      s=vprintf(s,'%*s: %4.1f*(MzzSa*Ia+MzzSy*Iy) - %4.1f*(MrrSa*Ia+MrrSy*Iy) = 0  1/w=%2.2e\n',ind,'CN',sin(LX.gpca(iS)),cos(LX.gpca(iS)),iWC(iS));
    end
    
    % Field gradient constraint
    if ~any(isnan([iWV(iS),LX.gpvrr(iS)]))
      if ~nx, ind = ind1; nx=true; else, ind=ind0; end
      s=vprintf(s,'%*s: BrrSa*Ia           = %5.2g                           1/w=%2.2e\n',ind,'dr B0r',LX.gpvrr(iS),iWV(iS));
    end
    if ~any(isnan([iWV(iS),LX.gpvrz(iS)]))
      if ~nx, ind = ind1; nx=true; else, ind=ind0; end
      s=vprintf(s,'%*s: BrzSa*Ia           = %5.2g                           1/w=%2.2e\n',ind,'dz B0r',LX.gpvrz(iS),iWV(iS));
    end
    if ~any(isnan([iWV(iS),LX.gpvzz(iS)]))
      if ~nx, ind = ind1; nx=true; else, ind=ind0; end
      s=vprintf(s,'%*s: BzzSa*Ia           = %5.2g                           1/w=%2.2e\n',ind,'dz Bz0',LX.gpvzz(iS),iWV(iS));
    end
  end
  if isempty(iS), s=vprintf(s,'   none\n'); end
end

if nargout==1
  varargout{1}=s;
end

end

function s=vprintf(s,varargin)
if iscell(s)
  s{end+1} = sprintf(varargin{:});
else
  fprintf(varargin{:})
end
end