/* [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved. */
# include "meq.h"
# ifdef SINGLE
# define FCT  sfbnd
# else
# define FCT  dfbnd
# endif

# define MAX(N,F,FB,K) for (i = N, pF = F; i--; pF++) if (*pF > FB) { FB = *pF; K = i; }
# define MIN(N,F,FB,K) for (i = N, pF = F; i--; pF++) if (*pF < FB) { FB = *pF; K = i; }

void FCT(FLT *FB, FLT *rB, FLT *zB, bool *lB, bool *lX,
         FLT *Fl, FLT *rl, FLT *zl,
         FLT *FX, FLT *rX, FLT *zX, int nX, FLT FN, int nl) {
 FLT FBX = FN, FBl = FN, *pF;
 int i, kX = -1, kl = -1;
 if (FN > 0.0) { /* Ip < 0 */
  MIN(nX, FX, FBX, kX)
  MIN(nl, Fl, FBl, kl)
  *lX = FBX < FBl;
 } else { /* Ip > 0 */
  MAX(nX, FX, FBX, kX)
  MAX(nl, Fl, FBl, kl)
  *lX = FBX > FBl;
 }
 if (++kX && *lX) {
  kX = nX - kX;
  *FB = FBX;
  *rB = rX[kX];
  *zB = zX[kX];
  *lB = 1;
 } else if (++kl) {
  kl = nl - kl;
  *FB = FBl;
  *rB = rl[kl];
  *zB = zl[kl];
  *lB = 1;
 } else {
  *FB = FN;
  *rB = *zB = FLTC(0.0);
  *lB = 0;
 }
}
