/* [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved. */
# include "meq.h"

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

 /* l = bavxmex(rx,zx,vrx,vzx,a,r,z,mode) */

# define MEXNAME bavxmex

# define L    pout[0]

# define RX   prhs[0]  /* double, def size=nx */
# define ZX   prhs[1]  /* double, size=nx     */
# define VRX  prhs[2]  /* double, size=nx     */
# define VZX  prhs[3]  /* double, size=nx     */
# define A    prhs[4]  /* numeric, scalar or empty */
# define R    prhs[5]  /* double, def size=[nr,nz] OR def size=nr */
# define Z    prhs[6]  /* double, size=nr*nz       OR def size=nz */

 CHECK_NARGIN_GE(7);
 CHECK_NARGIN_LE(8);
 
 int nX = mxGetNumberOfElements(RX);

 CHECK_DOUBLE(RX);
 CHECK_DOUBLE(ZX);
 CHECK_NUMEL(ZX,nX);
 CHECK_DOUBLE(VRX);
 CHECK_NUMEL(VRX,nX);
 CHECK_DOUBLE(VZX);
 CHECK_NUMEL(VZX,nX);
 CHECK_NUMERIC(A);
 CHECK_SCALAR_OR_EMPTY(A);
 CHECK_DOUBLE(R);
 CHECK_DOUBLE(Z);

 CHECK_NARGOUT_LE(1);

 mxArray *pout[1];

 double a = mxIsEmpty(A) ? 0.0 : mxGetScalar(A);
 if (nrhs == 7) {
  int nr = mxGetN(R);
  int nz = mxGetM(R);
  CHECK_NUMEL_MSG(Z,nr*nz," when called with 7 arguments");
  bool *l = (bool *)mxGetLogicals(L = mxCreateLogicalMatrix(nz,nr));
  dbavx(l,a,mxGetPr(RX),mxGetPr(ZX),mxGetPr(VRX),mxGetPr(VZX),nX,mxGetPr(R),mxGetPr(Z),nz*nr);
 } else {
  int nr = mxGetM(R)*mxGetN(R);
  int nz = mxGetM(Z)*mxGetN(Z);
  bool *l = (bool *)mxGetLogicals(L = mxCreateLogicalMatrix(nz,nr));
  dbavx2(l,a,mxGetPr(RX),mxGetPr(ZX),mxGetPr(VRX),mxGetPr(VZX),nX,mxGetPr(R),mxGetPr(Z),nr,nz);
 }

 ASSIGN_PLHS;
}
