/* [+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            )
    [y1,y2,y3,y4         ] = bf*mex(#,par, u1,F0,F1, q1, q2, q3) */

# define MEXNAME bfabmex

# 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 MXFA  prhs[3]  /* ... for remaining ... */
# define MXFB  prhs[4]  /* ... arguments in  ... */
# define q1    prhs[5]  /* ... separate mode ... */
# define q2    prhs[6]  /* ... blocks        ... */
# define q3    prhs[7]  /* ...               ... */

 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 = mxGetPr(par)[0], nT = mxGetPr(par)[1]; 
 if (nrhs>3)
 {
   CHECK_NUMERIC(MXFA);
   CHECK_NUMERIC(MXFB);
   if (!mxIsScalar(MXFA) || !mxIsScalar(MXFB)>1) {
     mexErrMsgIdAndTxt("bfabmex: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(nP+nT,1,mxREAL));
   double *fTg = mxGetPr(y2 = mxCreateDoubleMatrix(nP+nT,1,mxREAL));
   double *TDg = mxGetPr(y3 = mxCreateDoubleMatrix(1,nP+nT,mxREAL));
   switch (nP) {
    case 3: fPg[2] = 1.0; TDg[2] = 1.0;
    case 2: fPg[1] = 1.0; TDg[1] = 1.0;
    case 1: fPg[0] = 1.0; TDg[0] = 1.0;
   }
   switch (nT) {
    case 3: fTg[2+nP] = 1.0; TDg[2+nP] = 1.0;
    case 2: fTg[1+nP] = 1.0; TDg[1+nP] = 1.0;
    case 1: fTg[  nP] = 1.0; TDg[  nP] = 1.0;
   }
   break;
  }
  case 1 : {
   /* u1    prhs[2]  double, size=(nzy+2)*(nry+2) */
   /* MXFA  prhs[3]  numeric, scalar [checked] */
   /* MXFB  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);

   int ny = nzy*nry;
   double  *Tyg = mxGetPr(y1 = mxCreateDoubleMatrix(ny,nP+nT,mxREAL));
   double  *Tpg = mxGetPr(y2 = mxCreateDoubleMatrix(1 ,nP+nT,mxREAL));
   double *ITpg = mxGetPr(y3 = mxCreateDoubleMatrix(1 ,nP+nT,mxREAL));
   double  *Fx  =               mxGetPr(u1)+nzy+3;
   double   FB  =           mxGetScalar(MXFB);
   double   FAB =           mxGetScalar(MXFA) - FB;
   int8_t  *Opy =   (int8_t *)mxGetData(q1);
   double  *ry  =               mxGetPr(q2),
           *iry =               mxGetPr(q3);
   double  *TyP1 =  Tyg        ,  *TyP2 =  Tyg +        ny,  *TyP3 =  Tyg +  2    *ny,
           *TyT1 =  Tyg + nP*ny,  *TyT2 =  Tyg + (1+nP)*ny,  *TyT3 =  Tyg + (2+nP)*ny,
           *TpP1 =  Tpg        ,  *TpP2 =  Tpg +  1       ,  *TpP3 =  Tpg +  2       ,
           *TpT1 =  Tpg + nP   ,  *TpT2 =  Tpg +  1+nP    ,  *TpT3 =  Tpg +  2+nP    ,
          *ITpP1 = ITpg        , *ITpP2 = ITpg +  1       , *ITpP3 = ITpg +  2       ,
          *ITpT1 = ITpg + nP   , *ITpT2 = ITpg +  1+nP    , *ITpT3 = ITpg +  2+nP    ,
          g1,g2,g3,Ig1,Ig2,Ig3,g;
   int i,j;
   for (i=nry; i--; ry++, iry++, Fx+=2)
    for (j=nzy; j--; Fx++, TyP1++, TyP2++, TyP3++, TyT1++, TyT2++, TyT3++)
     if (*Opy++) {
       g1  = *Fx - FB ;
       g3  = g1 - FAB; /* Fx-FA */
       g2  = g1 * g3;
       g3  = g2 * (g1 + g3);
      Ig3  = g1 * g1;
      Ig1  = 0.5 * Ig3;
      Ig2  = Ig1 * (0.66666666666666666666 * g1 -FAB);
      Ig3 *= 0.5 * g2 * g2;
      switch (nP) {
       case 3:  *TyP3  = g =  *ry *  g3;
               *ITpP3 +=      *ry * Ig3;
                *TpP3 += g             ;
       case 2:  *TyP2  = g =  *ry *  g2;
               *ITpP2 +=      *ry * Ig2;
                *TpP2 += g             ;
       case 1:  *TyP1  = g =  *ry *  g1;
               *ITpP1 +=      *ry * Ig1;
                *TpP1 += g             ;
      }
      switch (nT) {
       case 3:  *TyT3  = g = *iry *  g3;
               *ITpT3 +=     *iry * Ig3;
                *TpT3 += g             ;
       case 2:  *TyT2  = g = *iry *  g2;
               *ITpT2 +=     *iry * Ig2;
                *TpT2 += g             ;
       case 1:  *TyT1  = g = *iry *  g1;
               *ITpT1 +=     *iry * Ig1;
                *TpT1 += g             ;
      }
     }
   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);
   double  *gQg = mxGetPr(y1 = mxCreateDoubleMatrix(nQ,nP+nT,mxREAL));
   double *IgQg = mxGetPr(y2 = mxCreateDoubleMatrix(nQ,nP+nT,mxREAL));
   double  *FNQ  =               mxGetPr(u1);
   double  *gQP1 =  gQg        ,  *gQP2 =  gQg +        nQ,  *gQP3 =  gQg +  2    *nQ,
           *gQT1 =  gQg + nP*nQ,  *gQT2 =  gQg + (1+nP)*nQ,  *gQT3 =  gQg + (2+nP)*nQ,
          *IgQP1 = IgQg        , *IgQP2 = IgQg +        nQ, *IgQP3 = IgQg +  2    *nQ,
          *IgQT1 = IgQg + nP*nQ, *IgQT2 = IgQg + (1+nP)*nQ, *IgQT3 = IgQg + (2+nP)*nQ,
          g1,g2,g3,Ig1,Ig2,Ig3,g;
   int i;
   for (i=nQ; i--;) {
     g2  = *FNQ++;
     g1  = g2 - 1.0;
     g3  = g2 * 2.0;
    Ig1  = 0.5 * g1 * g1;
    Ig2  = Ig1 * 0.33333333333333333333 * (g3 + 1.0);
    Ig3  = Ig1 *g2 * g2;
     g2 *= g1;
     g3  = g2 * (g3 - 1.0);
     switch (nP) {
      case 3:  *gQP3++ =  g3;
              *IgQP3++ = Ig3;
      case 2:  *gQP2++ =  g2;
              *IgQP2++ = Ig2;
      case 1:  *gQP1++ =  g1;
              *IgQP1++ = Ig1;
     }
     switch (nT) {
      case 3:  *gQT3++ =  g3;
              *IgQT3++ = Ig3;
      case 2:  *gQT2++ =  g2;
              *IgQT2++ = Ig2;
      case 1:  *gQT1++ =  g1;
              *IgQT1++ = Ig1;
     }
   }
   break;
  }
  case 3: {
   /* u1    prhs[2]  double, size=nP+nT */
   /* MXFA  prhs[3]  numeric, scalar [checked] */
   /* MXFB  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);

   double  *aPpg = mxGetPr(y1 = mxCreateDoubleMatrix(nP+nT,1,mxREAL));
   double *aTTpg = mxGetPr(y2 = mxCreateDoubleMatrix(nP+nT,1,mxREAL));
   double   *aPg = mxGetPr(y3 = mxCreateDoubleMatrix(nP+nT,1,mxREAL));
   double *ahqTg = mxGetPr(y4 = mxCreateDoubleMatrix(nP+nT,1,mxREAL));
   double  *ag   = mxGetPr(u1);
   double   FBA1 = mxGetScalar(MXFB) - mxGetScalar(MXFA),
            FBA2 = FBA1*FBA1, FBA3 = FBA2*FBA1;
   double  ids   = mxGetScalar(q3);
   double  cP    = ids*0.159154943091895 /* 1/2/pi */,
           cT    = ids*2e-7 /* mu0/2/pi */;
   switch (nP) {
    case 3:   aPg[2   ] = (  aPpg[2   ] = ag[2   ] * cP * FBA3 ) * FBA1;
    case 2:   aPg[1   ] = (  aPpg[1   ] = ag[1   ] * cP * FBA2 ) * FBA1;
    case 1:   aPg[0   ] = (  aPpg[0   ] = ag[0   ] * cP * FBA1 ) * FBA1;
   }
   switch (nT) {
    case 3: ahqTg[2+nP] = ( aTTpg[2+nP] = ag[2+nP] * cT * FBA3 ) * FBA1;
    case 2: ahqTg[1+nP] = ( aTTpg[1+nP] = ag[1+nP] * cT * FBA2 ) * FBA1;
    case 1: ahqTg[  nP] = ( aTTpg[  nP] = ag[  nP] * cT * FBA1 ) * FBA1;
   }
   break;
  }
  case 5: {
   /* u1    prhs[2]  unused */
   /* MXFA  prhs[3]  numeric, scalar [checked] */
   /* MXFB  prhs[4]  numeric, scalar [checked] */
   CHECK_NARGIN_EQ(5);
   CHECK_NARGOUT_LE(2);

   double  *gA = mxGetPr(y1 = mxCreateDoubleMatrix(1,nP+nT,mxREAL));
   double *IgA = mxGetPr(y2 = mxCreateDoubleMatrix(1,nP+nT,mxREAL));
   double FAB1 = mxGetScalar(MXFA) - mxGetScalar(MXFB),
          FAB2 = FAB1*FAB1, FAB3 = FAB2*FAB1;
   switch (nP) {
    case 3:  
    case 2:                IgA[   1] = -0.166666666666666*FAB3;
    case 1: gA[0 ] = FAB1; IgA[   0] =  0.500000000000000*FAB2;
   }
   switch (nT) {
    case 3:
    case 2:                IgA[1+nP] = -0.166666666666666*FAB3;
    case 1: gA[nP] = FAB1; IgA[  nP] =  0.500000000000000*FAB2;
   }
   break;
  }
  case 6: {
   /* u1    prhs[2]  unused */
   /* MXFA  prhs[3]  numeric, scalar [checked] */
   /* MXFB  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 ? nP-1 : 0, nTq = nT ? nT-1 : 0, nq = nPq + nTq;
   double *Qqg = mxGetPr(y1 = mxCreateDoubleMatrix(nq,nP+nT,mxREAL));
                         y2 = mxCreateDoubleMatrix(nq,    1,mxREAL) ;
   double FBA = mxGetScalar(MXFB)-mxGetScalar(MXFA), FBA2 = FBA*FBA, FBA3 = FBA2*FBA;
   double rA = mxGetScalar(q1), irA = mxGetScalar(q2), ids = mxGetScalar(q3);
   switch (nP) { case 3: Qqg[     2*nq+    1] = 3.464101615137754*FBA3*rA*ids; /* sqrt(12) */
                 case 2: Qqg[       nq      ] = 2.000000000000000*FBA2*rA*ids;
   }
   switch (nT) { case 3: Qqg[(nP+2)*nq+nPq+1] = 3.464101615137754*FBA3*irA*ids;
                 case 2: Qqg[(nP+1)*nq+nPq  ] = 2.000000000000000*FBA2*irA*ids;
   }   
   break;
  }
  case 7: {
   /* u1    prhs[2]  unused */
   /* MXFA  prhs[3]  numeric, scalar [checked] */
   /* MXFB  prhs[4]  numeric, scalar [checked] */
   CHECK_NARGIN_EQ(5);
   CHECK_NARGOUT_LE(2);

   int nPc = nP==3 ? 2 : nP, nTc = 0, nc = nPc + nTc;
   double *Qcg = mxGetPr(y1 = mxCreateDoubleMatrix(nc,nP+nT,mxREAL));
                         y2 = mxCreateDoubleMatrix(nc,    1,mxREAL) ;
   double FBA = mxGetScalar(MXFB)-mxGetScalar(MXFA);
   switch (nP) { case 3: Qcg[5] = FBA*FBA;
                 case 2: Qcg[2] = Qcg[3] = FBA; Qcg[1] = 1.0;
                 case 1: Qcg[0] = 3.0;
   }
   break;
  }
 }

 ASSIGN_PLHS;
}
