function [rA,zA,FA,dr2FA,dz2FA,drzFA,rX,zX,FX,dr2FX,dz2FX,drzFX,rl,zl,...
  rB,zB,FB,lB,lX,Opy] = meqpdom(Fx,Ip,isaddl,L)
%MEQPDOM  Plasma domain parameters
% [RA,ZA,FA,DR2FA,DZ2FA,DRZFA,RX,ZX,FX,DR2FX,DZ2FX,DRZFX,RL,ZL,
% RB,ZB,FB,LB,LX,OPY] = MEQPDOM(FX,IP,ISADDL,L) 
%
% MEQPDOM returns position, flux gradients of magnetic axes, x points,
% boundary points, as well as plasma domain parameters.
%
% FX(Z,R) is the flux map 
% IP plasma current for its sign
% ISADDL selects the use of X-points 0: exclude x-points 1: use x-points
%
% RA,ZA,FA position and flux of the magnetic axes
% DR2FA,DZ2FA,DRZFA flux 2nd derivatives on axes, 
% RX,ZX,FX position and flux of the X-points on the X-points polygon,
% DR2FX,DZ2FX,DRZFX flux 2nd derivatives at X-point, 
% RL,ZL the limiter points outside the divertor private zone, 
% RB,ZB,FB point and flux defining LCFS or internal separatrices
% LB is FALSE if no valid point defines the LCFS 
% LX is TRUE if an X-point defines the LCFS or separatrix
% OPY(Z,R) indices indicate different plasma domain.
%
% For details, see: [J-M. Moret et al. Fus.Eng.Des 2015], Section 2.3
% Multi-domain case for doublets extension added later
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

sIp = sign(Ip);
Opy = zeros(L.nzy,L.nry,'int8');
rl = L.G.rl; zl = L.G.zl;
rx = L.G.rx; zx = L.G.zx;

FN = L.FN*sIp;

% find all x/o points
[zA,rA,FA,dz2FA,dr2FA,drzFA,zX,rX,FX,dz2FX,dr2FX,drzFX,stat] = ...
  asxymex(Fx,L.G.zx,L.G.rx,L.P.dasm,L.dzx,L.drx,L.idzx,L.idrx,L.Oasx,L.dimw{:});

if (~stat)
  % Issue warning only, errors handled if nA~=1.
  % Note: should never happen if itert=0 since dimw=numel(L.Oly)
  s=warning('backtrace','off');
  warning('asxymex:workingSpace','Working space (%s) too small, increase DIMW',ifc(isempty(L.dimw),'[]',num2str(L.dimw{1})));
  warning(s);
end

if L.P.icsint
  [zA,rA,FA,dz2FA,dr2FA,drzFA,zX,rX,FX,dz2FX,dr2FX,drzFX] = ...
    asxycs(Fx,zA,rA,zX,rX,L);
end

nA = int32(numel(FA)); nX = int32(numel(FX));
switch nA
  case 0
    [rl,zl]    = deal([]   );
    [rB,zB,FB] = deal([]  );
    [lB,lX]    = deal(false);
    Opy = zeros(L.nzy,L.nry,'int8');
    return
  case 1
    vrX = rA - rX;
    vzX = zA - zX;
    %% LCFS flux
    if L.P.ilim == 0
      Fl = repmat(FN,L.G.nl,1);
      rl = L.G.rl; zl = L.G.zl;
    elseif L.P.ilim == 1 % default limiter treatment - only on rl,zl points
      Fl = fl4pmex(Fx,int32(L.kxl-1),L.clx,FN,L.G.rl,L.G.zl,rA,zA); % interpolate and keep only compatible sign of extrema
      Fl(~bavxmex(rX,zX,vrX,vzX,L.P.xdoma,L.G.rl,L.G.zl)) = FN; % exclude points outside x-point polygon
      rl = L.G.rl; zl = L.G.zl;
    elseif L.P.ilim == 2 % interpolated flux extrema on limiter
      iFl = bavxmex(rX,zX,vrX,vzX,L.P.xdoma,L.G.rl,L.G.zl); % include only points inside x-point polygon
      % interpolate on linear segments and keep only compatible sign of extrema
      [rl,zl,Fl] = fl4pinterp(Fx,L.kxl,L.clx,L.kxlh,L.clhx,iFl,FN,L.G.rl,L.G.zl,rA,zA);
    elseif L.P.ilim == 3 % Use cubic spline for Fx interpolation and pp interpolation for the limiter
      [Fl,rl,zl] = flcs(Fx,FN,rA,zA,L);
      Fl(~bavxmex(rX,zX,vrX,vzX,L.P.xdoma,rl,zl)) = FN; % exclude points outside x-point polygon
    else
      error('unknown ilim value %s',L.P.ilim)
    end

    if ~isaddl % no X points
      k = [];
    elseif isaddl==1 % keep X points on X points polygon
      k = bavxmex(rX,zX,vrX,vzX,0,rX,zX);
    else
      error('invalid value for isaddl %d',isaddl)
    end
    [FB,rB,zB,lB,lX] = fbndmex(Fl,rl,zl,FX(k),rX(k),zX(k),FN);
    
    % output valid limiter candidates
    rl = rl(Fl ~= FN);
    zl = zl(Fl ~= FN);
    %% Plasma domain
    Opy = int8(pdommex(Fx,FB,Ip,L.Oly & bavxmex(rX,zX,vrX,vzX,L.P.xdoma,L.ry,L.zy,[])));
    return
  otherwise
    % Multiple axes, multiple domains
    
    if ~L.P.ihole
      % Eliminate axes with wrong hessian, assuming all axes should have similar signs for current.
      iAok = (-sign(dr2FA+dz2FA) == sIp);
      FA=FA(iAok); rA=rA(iAok); zA=zA(iAok); dr2FA=dr2FA(iAok); dz2FA=dz2FA(iAok); drzFA=drzFA(iAok);
      nA = int32(numel(FA));
    end
    
    % order FA by vertical position, FX by flux value
    [zA,k] = sort(zA,'descend'); FA = FA(k); rA = rA(k); dr2FA = dr2FA(k); dz2FA = dz2FA(k); drzFA = drzFA(k);
    [FX,k] = sort(FX,'descend'); zX = zX(k); rX = rX(k); dr2FX = dr2FX(k); dz2FX = dz2FX(k); drzFX = drzFX(k);
    
    % Use X point hessians to compute directions of X point
    % (towards positive gradient direction)
    H1 = sIp*(dr2FX - dz2FX) * 0.5;
    H0 = sqrt( H1.*H1 + drzFX.*drzFX);
    vrX0 = sqrt(H0 + H1);
    vzX0 = sqrt(H0 - H1) .* ((drzFX >= 0)*2 -1)*sIp;
        
    %% Determine domain for each of the regions
    kO = (1:nX)';
    [kX,FB,rB,zB,FI,rI,zI,lX] = deal(zeros(1,nA+nX));
    lX = logical(lX); lB = false; np=[]; % init
    
    for kp = 1:nA+nX
      if kp <= nA
        % First treat domains that contain an extremum (axis)
        rI(kp) = rA(kp); zI(kp) = zA(kp); FI(kp) = FA(kp);
        
        vrX = rI(kp)-rX;
        vzX = zI(kp)-zX;
      else
        % Select first X-point connect to 2 already identified domains
        l = find(sum(kO == kX,2)==2,1);
        if isempty(l), break, end % If none is found, all remaining X-points are connected to a single domain and their separatix will intersect the wall. So our job is done.
        k = kO(l); kO(l) = []; % Eliminate selected X-point from remaining candidates.
        % Select an X point as flux limiting point for a
        % domain that does not contain a magnetic axis.
        rI(kp) = rX(k); zI(kp) = zX(k); FI(kp) = FX(k);
        vrX = vrX0;
        vzX = vzX0;
      end
      
      % use bavxmex to 'cut' regions outside of X-point polygon
      vsX = sign( vrX.*(rI(kp)-rX) + vzX.*(zI(kp)-zX));
      vrX = vrX .* vsX;
      vzX = vzX .* vsX;
      % cut regions to be excluded on (r,z)
      OX = bavxmex(rX(kO),zX(kO),vrX(kO),vzX(kO),0,rx,zx,[]);
      k = bavxmex(rX(kO),zX(kO),vrX(kO),vzX(kO),0,rX(kO),zX(kO)); % Eliminating X-point candidates outside X-point polygon
      
      %% LCFS flux
      if L.P.ilim == 0
        Fl = repmat(FN,L.G.nl,1);
        rl = L.G.rl; zl = L.G.zl;
      elseif L.P.ilim == 1 % default limiter treatment - only on rl,zl points
        Fl = fl4pmex(Fx,int32(L.kxl-1),L.clx,FN,L.G.rl,L.G.zl,rI(kp),zI(kp)); % interpolate and keep only compatible sign of extrema
        Fl(~bavxmex(rX,zX,vrX,vzX,L.P.xdoma,L.G.rl,L.G.zl)) = FN; % exclude points outside x-point polygon
        rl = L.G.rl; zl = L.G.zl;
      elseif L.P.ilim == 2 % interpolated flux extrema on limiter
        iFl = bavxmex(rX,zX,vrX,vzX,L.P.xdoma,L.G.rl,L.G.zl); % include only points inside x-point polygon
        % interpolate on linear segments and keep only compatible sign of extrema
        [rl,zl,Fl] = fl4pinterp(Fx,L.kxl,L.clx,L.kxlh,L.clhx,iFl,FN,L.G.rl,L.G.zl,rI(kp),zI(kp));
      elseif L.P.ilim == 3
        [Fl,rl,zl] = flcs(Fx,FN,rI(kp),zI(kp),L);
        Fl(~bavxmex(rX,zX,vrX,vzX,L.P.xdoma,rl,zl)) = FN; % exclude points outside x-point polygon
      else
        error('unknown ilim value %s',L.P.ilim)
      end
      
      % Candidate limiter points (flux extrema on limiter)
      [FB(kp),rB(kp),zB(kp),lB,lX(kp)] = fbndmex(Fl,rl,zl,FX(kO(k)),rX(kO(k)),zX(kO(k)),FN);
      rl = rl(Fl~=FN); zl = zl(Fl~=FN);  % output only valid limiter points (admittedly for last kp iteration only)

      if ~lB, break; end
      
      if lX(kp), kX(kp) = find(FX == FB(kp)); end % If domain is X-point limited, add its index to the list of boundary X-points.
      
      Oxyk = OX & (Fx -FI(kp)).*(Fx - FB(kp)) < 0;

      Opy(Oxyk(L.lxy)) = kp; % indicized Opy
      
      np = kp; % % update number of domains
    end
    % truncate to effective number of domains
    FB = FB(1:np)'; rB = rB(1:np)'; zB = zB(1:np)'; lX = lX(1:np)';
   
    % mask outside limiter
    Opy(~L.Oly) = 0;
end
end
