/* [+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[])
{
 
 /* [Opy] = pdommex(Fx,FB,Ip,OXy) */

# define MEXNAME pdommex

# define OPY  pout[0]

# define FX   prhs[0] /* double, def size=[nzy+2,nry+2] */
# define FB   prhs[1] /* numeric, scalar */
# define IP   prhs[2] /* numeric, scalar */
# define OXY  prhs[3] /* logical, size=[nzy,nry] */

 CHECK_NARGIN_EQ(4);

 int nzy = mxGetM(FX)-2, nry = mxGetN(FX)-2;

 CHECK_DOUBLE(FX);
 CHECK_NUMERIC(FB);
 CHECK_SCALAR(FB);
 CHECK_NUMERIC(IP);
 CHECK_SCALAR(IP);
 CHECK_LOGICAL(OXY);
 CHECK_NUMEL(OXY,nzy*nry);

 CHECK_NARGOUT_LE(1);

 mxArray *pout[1];

 OPY = mxCreateLogicalMatrix(nzy,nry);
 
 dpdom((bool *)mxGetLogicals(OPY), mxGetPr(FX), mxGetScalar(FB), mxGetScalar(IP), (bool *)mxGetLogicals(OXY), nzy, nry);

 ASSIGN_PLHS;
}
