/* [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved. */
# include "meq.h"
# ifdef SINGLE
# define FCT  sfl4p
# define BINT sbint
# else
# define FCT  dfl4p
# define BINT dbint
# endif
# define TEST if (df1*df2 > FLTC(0.0) | (df1-df2)*fn > FLTC(0.0)) *fl = fn;

void FCT(FLT *fl, FLT *f, int *k, FLT *c, FLT fn, FLT *rl, FLT *zl, FLT ra, FLT za, int nl, int nz) {
 int l;
 FLT f0, df1, df2, *fl0 = fl, *c0 = c;
 int *k0 = k;
 FLT drfl, dzfl, dafl;
 /* Interpolate flux on limiter points */
 BINT(fl, f, k, c, nl, nz);
 if (fn != FLTC(0.0)) {
  /* Discard non-extremum points */
  fl = fl0;
  f0 = *fl;
  df1 = f0 - fl[nl-1];
  df2 = fl[1] - f0;
  for (l=nl-2; l--; ) {
   TEST; fl++;
   df1 = df2;
   df2 = fl[1] - *fl;
  }
  TEST; fl++;
  df1 = df2;
  df2 = f0 - *fl;
  TEST;
  /* Discard points with incorrect sign of derivative
     in the direction away from plasma axis */
  fl = fl0; /* Restore pointers to     */
  k  = k0;  /*                         */
  c  = c0;  /* their initial locations */
  for (l=0; l < nl; l++, k++, c+=4, fl++) {
   if (*fl != fn) {
    drfl = (f[(*k)+nz  ] - f[(*k)     ])*(c[1]+c[0]) +
           (f[(*k)+nz+1] - f[(*k)+1   ])*(c[3]+c[2]);
    dzfl = (f[(*k)+1   ] - f[(*k)     ])*(c[3]+c[0]) +
           (f[(*k)+nz+1] - f[(*k)+nz  ])*(c[2]+c[1]);
    dafl = drfl*(ra-rl[l])+dzfl*(za-zl[l]);
    if (dafl*fn > 0) *fl = fn;
   }
  }
 }
}
