/* [+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[])
{
  
 /* [aS] = locSmex(aQ, qQ, qS, NF); */

# define MEXNAME locSmex

# define AS   pout[0]

# define AQ   prhs[0]  /* double, def size=[nO,nQ] */
# define QQ   prhs[1]  /* double, size=nQ */
# define QS   prhs[2]  /* double, def size=nS */
# define NF   prhs[3]  /* numeric, scalar, optional */

 CHECK_NARGIN_GE(3);
 CHECK_NARGIN_LE(4);

 int nO = mxGetM(AQ), nQ = mxGetN(AQ), nS = mxGetNumberOfElements(QS);

 CHECK_DOUBLE(AQ);
 CHECK_DOUBLE(QQ);
 CHECK_NUMEL(QQ,nQ);
 CHECK_DOUBLE(QS);
 if (nrhs == 4) { CHECK_NUMERIC(NF);
                  CHECK_SCALAR(NF); }

 CHECK_NARGOUT_LE(1);

 mxArray *pout[1];

 double nf;

 if (nrhs == 4) nf = mxGetScalar(NF);
 else           nf = NAN;
 
 AS = mxCreateDoubleMatrix(nO,nS,mxREAL);
 
 dlocQ(mxGetPr(AS), mxGetPr(AQ), mxGetPr(QQ), mxGetPr(QS), nf, nQ, nO, nS, 1);

 ASSIGN_PLHS;
}
