/* [+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[])
{

 /* [aR] = locRmex(aQ, qQ, qR[, naR[, NF]]); */

# define MEXNAME locRmex

# define AR   pout[0]

# define AQ   prhs[0] /* double, def size=nQ */
# define QQ   prhs[1] /* double, size=nQ */
# define QR   prhs[2] /* double, def size=nR */
# define NAR  prhs[3] /* numeric, scalar, optional */
# define NF   prhs[4] /* numeric, scalar, optional */

 int nQ, nR, naR;
 double nf;
 
 CHECK_NARGIN_GE(3);
 CHECK_NARGIN_LE(5);

 nQ = mxGetNumberOfElements(AQ);
 nR = mxGetNumberOfElements(QR);
 
 CHECK_DOUBLE(AQ);
 CHECK_DOUBLE(QQ);
 CHECK_NUMEL(QQ,nQ);
 CHECK_DOUBLE(QR);
 if (nrhs >= 4) { CHECK_NUMERIC(NAR);
                  CHECK_SCALAR(NAR); }
 if (nrhs >= 5) { CHECK_NUMERIC(NF);
                  CHECK_SCALAR(NF); }

 CHECK_NARGOUT_LE(1);

 mxArray *pout[1];
 
 if (nrhs >= 5) nf = mxGetScalar(NF);
 else           nf = NAN;

 if (nrhs >= 4) naR = mxGetScalar(NAR);
 else           naR = 1;
 
 AR = mxCreateDoubleMatrix(nR,naR,mxREAL);
 
 dlocQ(mxGetPr(AR), mxGetPr(AQ), mxGetPr(QQ), mxGetPr(QR), nf, nQ, 1, nR, naR);

 ASSIGN_PLHS;
}
