/* [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved. */
# include "meq.h"
# ifdef SINGLE
# define ABS  fabsf
# define FCT  sbavx
# define FCT2 sbavx2
# else
# define ABS  fabs
# define FCT  dbavx
# define FCT2 dbavx2
# endif

void FCT(bool *l, FLT a,  FLT *rx, FLT *zx, FLT *vrx, FLT *vzx,int nx, FLT *r, FLT *z, int n) {
 FLT *rxp, *zxp, *vrxp, *vzxp, rv, zv, rpx, zpx, vr, vz;
 int kx, k;
 bool ll = false /* in case nx = 0 */;
 if (a == FLTC(0.0))
  for (k=n; k--; ) {
   rv = *r++;
   zv = *z++;
   for (kx=nx, rxp=rx, zxp=zx, vrxp=vrx, vzxp=vzx; kx--; ) {
    if ((ll = (*vrxp++ * (rv - *rxp++) + *vzxp++ * (zv - *zxp++)) < FLTC(0.0))) break;
   }
   *l++ = !ll;
  }
 else
  for (k=n; k--; ) {
   rv = *r++;
   zv = *z++;
   for (kx=nx, rxp=rx, zxp=zx, vrxp=vrx, vzxp=vzx; kx--; ) {
    rpx = rv - *rxp++;
    zpx = zv - *zxp++;
    vr = *vrxp++;
    vz = *vzxp++;
    if ((ll = (vr*rpx + vz*zpx) < a * ABS(vz*rpx - vr*zpx))) break;
   }
   *l++ = !ll;
  }
}

void FCT2(bool *l, FLT a, FLT *rx, FLT *zx, FLT *vrx, FLT *vzx, int nx, FLT *r, FLT *z, int nr, int nz) {
 FLT *rxp, *zxp, *vrxp, *vzxp, *zp, rv, zv, rpx, zpx, vr, vz;
 int kx, kr, kz;
 bool ll = false /* in case nx = 0 */;
 if (a == FLTC(0.0))
  for (kr=nr; kr--;) {
   rv = *r++;
   for (kz=nz, zp=z; kz--; ) {
    zv = *zp++;
    for (kx=nx, rxp=rx, zxp=zx, vrxp=vrx, vzxp=vzx; kx--; ) {
     if ((ll = (*vrxp++ * (rv - *rxp++) + *vzxp++ * (zv - *zxp++)) < FLTC(0.0))) break;
    }
    *l++ = !ll;
   }
  }
 else
  for (kr=nr; kr--;) {
   rv = *r++;
   for (kz=nz, zp=z; kz--; ) {
    zv = *zp++;
    for (kx=nx, rxp=rx, zxp=zx, vrxp=vrx, vzxp=vzx; kx--; ) {
     rpx = rv - *rxp++;
     zpx = zv - *zxp++;
     vr = *vrxp++;
     vz = *vzxp++;
     if ((ll = (vr*rpx + vz*zpx) < a * ABS(vz*rpx - vr*zpx))) break;
    }
    *l++ = !ll;
   }
  }
}
