%MEQFW - Determine flux values to be tracked by gap algorithm
% FW = meqFW(FB,FX,Ip,nFW)
% Returns first the boundary value, then X points ordered by flux
% NB: This may include internal x-points in doublets
%
% Inputs: 
%   FB flux at boundary
%   FX fluxes at all X points
%   lX  diverted or not
%   nX number of X points (since FX may be larger than that)
%   nFW number of flux values to be tracked.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function FW = meqFW(FB,FX,lX,nX,nFW)
%#codegen

% List of values to track
FW = zeros(1,nFW,class(FB)); % same type as FB, important for slx
if nFW == 0, return; end
[~,isFX] = sort(abs(FX(1:nX) - FB)); sFX = FX(isFX); % sort FX from closest to furtherest from FB
if lX
  nFWX = min(nX,nFW);           % number of FW points that are X points
  FW(1:nFWX) = sFX(1:nFWX);     % fill sorted X points
  FW((nFWX+1):nFW)=sFX(nFWX);   % fill remaining ones with last FW
elseif nX>0                     % first is limiter point, then x points
  FW(1) = FB;
  nFWX = min(nX,nFW-1);         % number of FW points that are X points
  FW(1+(1:nFWX)) = sFX(1:nFWX); % fill sorted X points
  FW((nFWX+2):nFW) = FW(1+nFWX);% fill rest
else
  FW(1:nFW) = FB;               % no X points, track limiter boundary everywhere
end
end
