function LX = fbtx(L,LX)
% function LX = fbtx(L,LX)
% Defaults and general things for fbtx
%
% LX fields in FBT:
%
% Fields that should come from fbtxtok:
%  .t        Time grid
%  .IpD      Plasma current per domain
%  .bpD      betap per domain
%  .qA       Axis q per domain
%  .rBt      R*Bt
%  .gp*      Various equilibrium constraint/cost function parameters documented in FBTGP and FBTHELP
%  .idoublet This time slice is a doublet
%  .bfp      (optional) time-varying basis function parameters when using bffbt
%
% Fields that should come from fbtp or can be passed as arguments (see help FBTP):
%  .niter
%  .tol 
%  .issym 
%  .capaj
% 
% Fields from meqp:
%  .isaddl
%  .idoublet
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

nt = numel(LX.t);

%% defaults

pfields = {'idoublet','isaddl','issym','tol','capaj','niter'};
for field=pfields
  myfield = field{1};
  if isfield(LX,myfield)
    % if LX already has them
    assert(size(LX.(myfield),2)==numel(LX.t),...
      'LX.%s has incorrect number of time slices, expect %d found %d',myfield,numel(LX.t),size(LX.(myfield),2))
    continue % just check size
  else
    if isfield(L.P,myfield)
      % copy value from L.P.myfield
      pvalue = L.P.(myfield);
      if size(pvalue,2) == 1,         LX.(myfield) = repmat(pvalue,1,nt);
      elseif size(pvalue,2) == nt,    LX.(myfield) = pvalue;
      elseif size(pvalue,2) == numel(L.P.t) % if still on P.t time base from fbtptcv
                                      LX.(myfield) = pvalue(:,iround(L.P.t,LX.t)); 
      else, error('invalid size for L.P.%s. Expected single or %d time slices, found %d',...
          myfield,nt,size(pvalue,2))
      end
    else
      error('LX.%s is not provided nor can a default be found in L.P',myfield)
    end
  end
end

%% bffbt special
if ~isfield(LX,'bfp')
  bfct = func2str(L.bfct);
  if strcmp(bfct,'bffbt') || (strcmp(bfct,'bfgenD') && strcmp(func2str(L.bfp{1}),'bffbt'))
    LX.bfp = L.bfp;
  end
end
  
%% multidomain special
LX = meqxdoublet(L,LX);

%% Remove non-D inputs to avoid confusion
for field={'bp','Ip','bt'}
  if isfield(LX,field{:})
    LX = rmfield(LX,field{:});
  end
end

%% Merge constraints that can be merged
if L.P.mergex
  for kt=1:nt
    LXt = meqxk(LX,kt);
    iok = ~any(isnan([LXt.gpr,LXt.gpz]),2);
    % grouping Br=0 and Bz=0 constraints for X-points.
    xx = find(iok & ((LXt.gpbr==0 & LXt.gpbd*LXt.gpbe==0) | (LXt.gpbz==0 & LXt.gpbd*LXt.gpbe==0)));
    [C,~,iC] = unique([LXt.gpr(xx),LXt.gpz(xx)],'rows');
    for ii = 1:size(C,1)% for each unique one
      ioccurrence = xx(iC==ii);
      if numel(ioccurrence)>1 && all(any([LXt.gpbr(ioccurrence),LXt.gpbz(ioccurrence)]==0)) % more than one occurrence and both Br=0 and Bz=0
        ikeep = ioccurrence(1    ); % keep first
        ielim = ioccurrence(2:end); % eliminate others
        LX.gpbr(ikeep,kt) = 0; LX.gpbz(ikeep,kt) = 0;   LX.gpbe(ikeep,kt)=0;
        LX.gpbr(ielim,kt)=NaN; LX.gpbz(ielim,kt) = NaN; LX.gpbe(ielim,kt)=NaN;
      end
    end
  end
end

%% Group gp* by same gpr,gpz points
list = strcat('gp',{'r','z','b','fa','fb','fe','br','bz','ba','be','cr','cz','ca','ce','vrr','vrz','vzz','ve'});
for kt = 1:nt
  [~,ia,ic] = unique([LX.gpr(:,kt),LX.gpz(:,kt)],'rows'); % ia(ic) is the index of the first occurence of each row
  [~,sorti] = sort(ia(ic));                               % This reorders rows to group them by value while keeping the order stable
  for field = list
    LX.(field{1})(:,kt) = LX.(field{1})(sorti,kt);
  end
end

%% Checks
assert(size(LX.gpia,1)==L.G.na,'gpia size does not match G.na')
assert(size(LX.gpie,1)==L.G.na,'gpie size does not match G.na')

%% meqx for general additions and further defaults
LX = meqx(L,LX);

%% Verbosity
if L.P.debug>1
  fbtxdisp(L,LX);
end

end
