/* [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved. */
# include "genlib.h"

void mexFunction(int nlhs,       mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{

# define MEXNAME ztrig

 /* y = ztrig(x,xl,xh); */

# define Y   pout[0] 

# define X   prhs[0] /* double */
# define XL  prhs[1] /* numeric, scalar */
# define XH  prhs[2] /* numeric, scalar */

 CHECK_NARGIN_EQ(3);

 CHECK_DOUBLE(X);
 CHECK_NUMERIC(XL);
 CHECK_SCALAR(XL);
 CHECK_NUMERIC(XH);
 CHECK_SCALAR(XH);

 CHECK_NARGOUT_LE(1);

 mxArray *pout[1];

 register int i,isl,ish;
 int mX,nX;
 double xl,xh,*x,*y;

 mX = mxGetM(X); nX = mxGetN(X);
 Y = mxCreateDoubleMatrix(mX,nX,mxREAL);
 if (mX * nX == 0) return;

 x = mxGetPr(X); y = mxGetPr(Y);
 xl = mxGetScalar(XL); xh = mxGetScalar(XH);
 y[0] = (double)(x[0] <= ((xl + xh) / 2));
 for (i = 1; i < mX*nX; i++) {
  isl = x[i] < xl; ish = x[i] > xh;
  y[i] = (double)(ish || (int)y[i-1] && !(isl || ish));
 }

 ASSIGN_PLHS;
}
