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

 /* [fPg,fTg,TDg         ] = bf*mex(0,par                              )
    [Tyg,TpDg,ITpDg      ] = bf*mex(1,par, Fx,F0,F1,Opy, ry,iry        )
    [gQDg,IgQDg          ] = bf*mex(2,par,FNQ,F0,F1                    )
    [aPpg,aTTpg,aPg,ahqTg] = bf*mex(3,par, ag,F0,F1,fPg,fTg,ids        )
    [gA,IgA              ] = bf*mex(5,par, ~ ,F0,F1                    )
    [Qqg,Xq              ] = bf*mex(6,par, ~ ,F0,F1, r0,ir0,ids        )
    [Qcg,Xc              ] = bf*mex(7,par, ~ ,F0,F1, r0,ir0,ids        )
    [Bty                 ] = bf*mex(8,par, Fx,F0,F1,Opy, ag,rBt,ids,iry)
    [y1,y2,y3,y4         ] = bf*mex(#,par, u1,F0,F1, q1, q2, q3, q4, q5) */

# define MEXNAME bf3pmex

# define y1    pout[0]
# define y2    pout[1]
# define y3    pout[2]
# define y4    pout[3]

# define MODE  prhs[0]  /* numeric, scalar */
# define par   prhs[1]  /* logical, scalar */
# define u1    prhs[2]  /* ... see checks     ...*/
# define FA    prhs[3]  /* ... for remaining ... */
# define FB    prhs[4]  /* ... arguments in  ... */
# define q1    prhs[5]  /* ... separate mode ... */
# define q2    prhs[6]  /* ... blocks        ... */
# define q3    prhs[7]  /* ...               ... */
# define q4    prhs[8]  /* ...               ... */
# define q5    prhs[9]  /* ...               ... */

 CHECK_NARGIN_GE(2);
 CHECK_NUMERIC(MODE);
 CHECK_SCALAR(MODE);
 CHECK_LOGICAL(par);
 CHECK_SCALAR(par);

 mxArray *pout[4];

 int mode = mxGetScalar(MODE);
 if (nrhs>3)
 {
   CHECK_NUMERIC(FA);
   CHECK_NUMERIC(FB);
   if (!mxIsScalar(FA) || !mxIsScalar(FB)>1) {
     mexErrMsgIdAndTxt("bf3pmex:multiD","To use bf3pmex with multiple domains, please use bfgenD");
   }
 }
 switch (mode) {
  case 0 : {
   CHECK_NARGIN_EQ(2);
   CHECK_NARGOUT_LE(3);

   double *fPg = mxGetPr(y1 = mxCreateDoubleMatrix(3,1,mxREAL));
   double *fTg = mxGetPr(y2 = mxCreateDoubleMatrix(3,1,mxREAL));
   double *TDg = mxGetPr(y3 = mxCreateDoubleMatrix(1,3,mxREAL));
   fPg[0] = 1.0; fPg[1] = fPg[2] = 0.0;
   fTg[0] = 0.0; fTg[1] = fTg[2] = 1.0;
   TDg[0] = TDg[1] = TDg[2] = 1.0;
   break;
  }
  case 1 : {
   /* u1    prhs[2]  double, size=(nzy+2)*(nry+2) */
   /* FA    prhs[3]  numeric, scalar [checked] */
   /* FB    prhs[4]  numeric, scalar [checked] */
   /* q1    prhs[5]  int8, def size=[nzy,nry] */
   /* q2    prhs[6]  double, size=nry */
   /* q3    prhs[7]  double, size=nry */
   CHECK_NARGIN_EQ(8);
   int nzy = mxGetM(q1), nry = mxGetN(q1);
   CHECK_DOUBLE(u1);
   CHECK_NUMEL(u1,(nzy+2)*(nry+2));
   CHECK_INT8(q1);
   CHECK_DOUBLE(q2);
   CHECK_NUMEL(q2,nry);
   CHECK_DOUBLE(q3);
   CHECK_NUMEL(q3,nry);
   CHECK_NARGOUT_LE(3);

   y1 = mxCreateDoubleMatrix(nzy*nry,3,mxREAL);
   y2 = mxCreateDoubleMatrix(1      ,3,mxREAL);
   y3 = mxCreateDoubleMatrix(1      ,3,mxREAL);
   dbf3p1(mxGetPr(y1), mxGetPr(y2), mxGetPr(y3),
          mxGetPr(u1), mxGetScalar(FA), mxGetScalar(FB),
          (int8_t *)mxGetData(q1), mxGetPr(q2), mxGetPr(q3),
          nzy, nry,
          mxIsLogicalScalarTrue(par), 0);
   break;
  }
  case 2: {
   /* u1    prhs[2]  double, def size=nQ */
   /* FA    prhs[3]  unused */
   /* FB    prhs[4]  unused */
   CHECK_NARGIN_EQ(5);
   CHECK_DOUBLE(u1);
   CHECK_NARGOUT_LE(2);

   int nQ = mxGetNumberOfElements(u1);
   y1 = mxCreateDoubleMatrix(nQ,3,mxREAL);
   y2 = mxCreateDoubleMatrix(nQ,3,mxREAL);
   dbf3p2(mxGetPr(y1), mxGetPr(y2), mxGetPr(u1), nQ);
   break;
  }
  case 3: {
   /* u1    prhs[2]  double, size=3 */
   /* FA    prhs[3]  numeric, scalar [checked] */
   /* FB    prhs[4]  numeric, scalar [checked] */
   /* q1    prhs[5]  unused */
   /* q2    prhs[6]  unused */
   /* q3    prhs[7]  numeric, scalar */
   CHECK_NARGIN_EQ(8);
   CHECK_DOUBLE(u1);
   CHECK_NUMEL(u1,3);
   CHECK_NUMERIC(q3);
   CHECK_SCALAR(q3);
   CHECK_NARGOUT_LE(4);

   y1 = mxCreateDoubleMatrix(3,1,mxREAL);
   y2 = mxCreateDoubleMatrix(3,1,mxREAL);
   y3 = mxCreateDoubleMatrix(3,1,mxREAL);
   y4 = mxCreateDoubleMatrix(3,1,mxREAL);
   dbf3p3(mxGetPr(y1), mxGetPr(y2), mxGetPr(y3), mxGetPr(y4),
          mxGetPr(u1), mxGetScalar(FA), mxGetScalar(FB),
          mxGetScalar(q3));
   break;
  }
  case 5: {
   /* u1    prhs[2]  unused */
   /* FA    prhs[3]  numeric, scalar [checked] */
   /* FB    prhs[4]  numeric, scalar [checked] */
   CHECK_NARGIN_EQ(5);
   CHECK_NARGOUT_LE(2);

   y1 = mxCreateDoubleMatrix(1,3,mxREAL);
   y2 = mxCreateDoubleMatrix(1,3,mxREAL);
   dbf3p5(mxGetPr(y1), mxGetPr(y2), mxGetScalar(FA), mxGetScalar(FB),
    mxIsLogicalScalarTrue(par));
   break;
  }
  case 6: {
   /* u1    prhs[2]  unused */
   /* FA    prhs[3]  numeric, scalar [checked] */
   /* FB    prhs[4]  numeric, scalar [checked] */
   /* q1    prhs[5]  unused */
   /* q2    prhs[6]  numeric, scalar */
   /* q3    prhs[7]  numeric, scalar */
   CHECK_NARGIN_EQ(8);
   CHECK_NUMERIC(q2);
   CHECK_SCALAR(q2);
   CHECK_NUMERIC(q3);
   CHECK_SCALAR(q3);
   CHECK_NARGOUT_LE(2);

   y1 = mxCreateDoubleMatrix(1,3,mxREAL);
   y2 = mxCreateDoubleMatrix(1,1,mxREAL);
   dbf3p6(mxGetPr(y1), mxGetPr(y2), mxGetScalar(FA), mxGetScalar(FB),
     mxGetScalar(q2),  mxGetScalar(q3));
   break;
  }
  case 7: {
   /* u1    prhs[2]  unused */
   /* FA    prhs[3]  numeric, scalar [checked], unused */
   /* FB    prhs[4]  numeric, scalar [checked], unused */
   CHECK_NARGIN_EQ(5);
   CHECK_NARGOUT_LE(2);

   y1 = mxCreateDoubleMatrix(0,3,mxREAL);
   y2 = mxCreateDoubleMatrix(0,1,mxREAL);
   break;
  }
  case 8: {
   /* u1    prhs[2]  double, size=(nzy+2)*(nry+2) */
   /* FA    prhs[3]  numeric, scalar [checked] */
   /* FB    prhs[4]  numeric, scalar [checked] */
   /* q1    prhs[5]  int8, def size=[nzy,nry] */
   /* q2    prhs[6]  double, size=3 */
   /* q3    prhs[7]  numeric, scalar */
   /* q4    prhs[8]  numeric, scalar */
   /* q5    prhs[9]  double, size=nry */
   CHECK_NARGIN_EQ(10);
   int nzy = mxGetM(q1), nry = mxGetN(q1);
   CHECK_DOUBLE(u1);
   CHECK_NUMEL(u1,(nzy+2)*(nry+2));
   CHECK_INT8(q1);
   CHECK_DOUBLE(q2);
   CHECK_NUMEL(q2,3);
   CHECK_NUMERIC(q3);
   CHECK_SCALAR(q3);
   CHECK_NUMERIC(q4);
   CHECK_SCALAR(q4);
   CHECK_DOUBLE(q5);
   CHECK_NUMEL(q5,nry);
   CHECK_NARGOUT_LE(1);

   y1 = mxCreateDoubleMatrix(nzy,nry,mxREAL);
   dbf3p8(mxGetPr(y1), mxGetPr(u1), mxGetScalar(FA), mxGetScalar(FB), 
    (int8_t *)mxGetData(q1), mxGetPr(q2), mxGetScalar(q3), 
    mxIsLogicalScalarTrue(par), mxGetScalar(q4), mxGetPr(q5), nzy, nry);
   break;
  }
 }

 ASSIGN_PLHS;
}
