/* [+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        )
    [g0g,Ig0g            ] = 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                    )
    [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 bfefmex

# 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]  /* double, size=2  */
# 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_DOUBLE(par);
 CHECK_NUMEL(par,2);

 mxArray *pout[4];

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

   double *fPg = mxGetPr(y1 = mxCreateDoubleMatrix(nP+nT,1,mxREAL));
   double *fTg = mxGetPr(y2 = mxCreateDoubleMatrix(nP+nT,1,mxREAL));
   double *TDg = mxGetPr(y3 = mxCreateDoubleMatrix(1,nP+nT,mxREAL));
   int i;
   for (i=0; i<nP; i++) {
     fPg[   i] = 1.0; fTg[   i] = 0.0; TDg[   i] = 1.0;
   }
   for (i=0; i<nT; i++) {
     fPg[nP+i] = 0.0; fTg[nP+i] = 1.0; TDg[nP+i] = 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=nzy */
   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,nP+nT,mxREAL);
   y2 = mxCreateDoubleMatrix(1      ,nP+nT,mxREAL);
   y3 = mxCreateDoubleMatrix(1      ,nP+nT,mxREAL);
   dbfef1(mxGetPr(y1), mxGetPr(y2), mxGetPr(y3),
          mxGetPr(u1), mxGetScalar(FA), mxGetScalar(FB),
          (int8_t *)mxGetData(q1), mxGetPr(q2), mxGetPr(q3),
          nzy, nry,
          nP, nT, 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,nP+nT,mxREAL);
   y2 = mxCreateDoubleMatrix(nQ,nP+nT,mxREAL);
   dbfef2(mxGetPr(y1), mxGetPr(y2), mxGetPr(u1), nQ, nP, nT);
   break;
  }
  case 3: {
   /* u1    prhs[2]  double, size=nP+nT */
   /* 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,nP+nT);
   CHECK_NUMERIC(q3);
   CHECK_SCALAR(q3);
   CHECK_NARGOUT_LE(4);
  
   y1 = mxCreateDoubleMatrix(nP+nT,1,mxREAL);
   y2 = mxCreateDoubleMatrix(nP+nT,1,mxREAL);
   y3 = mxCreateDoubleMatrix(nP+nT,1,mxREAL);
   y4 = mxCreateDoubleMatrix(nP+nT,1,mxREAL);
   dbfef3(mxGetPr(y1), mxGetPr(y2), mxGetPr(y3), mxGetPr(y4),
          mxGetPr(u1), mxGetScalar(FA), mxGetScalar(FB),
          mxGetScalar(q3), nP, nT);
   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,nP+nT,mxREAL);
   y2 = mxCreateDoubleMatrix(1,nP+nT,mxREAL);
   dbfef5(mxGetPr(y1), mxGetPr(y2), mxGetScalar(FA), mxGetScalar(FB),
          nP, nT);
   break;
  }
  case 6: {
   /* u1    prhs[2]  unused */
   /* FA    prhs[3]  numeric, scalar [checked] */
   /* FB    prhs[4]  numeric, scalar [checked] */
   /* q1    prhs[5]  numeric, scalar */
   /* q2    prhs[6]  numeric, scalar */
   /* q3    prhs[7]  numeric, scalar */
   CHECK_NARGIN_EQ(8);
   CHECK_NUMERIC(q1);
   CHECK_SCALAR(q1);
   CHECK_NUMERIC(q2);
   CHECK_SCALAR(q2);
   CHECK_NUMERIC(q3);
   CHECK_SCALAR(q3);
   CHECK_NARGOUT_LE(2);

   int nPq = nP>2 ? nP-2 : 0, nTq = nT>2 ? nT-2 : 0, nq = nPq + nTq + 2*((nP+nT)>0);
   y1 = mxCreateDoubleMatrix(nq,nP+nT,mxREAL);
   y2 = mxCreateDoubleMatrix(nq,1    ,mxREAL);
   dbfef6(mxGetPr(y1), mxGetPr(y2), mxGetScalar(FA), mxGetScalar(FB),
                mxGetScalar(q1), mxGetScalar(q2), mxGetScalar(q3), nP, nT, nq);
   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,nP+nT,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=nP+nT */
   /* 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,nP+nT);
   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);
   dbfef8(mxGetPr(y1), mxGetPr(u1), mxGetScalar(FA), mxGetScalar(FB), 
          (int8_t *)mxGetData(q1), mxGetPr(q2), mxGetScalar(q3), 
          nP, nT, mxGetScalar(q4), mxGetPr(q5), nzy, nry);
   break;
  }
 }

 ASSIGN_PLHS;
}
