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

 /* [rgeom,zgeom,amino,epsil,kappa,delta,deltl,deltu,rrmax,zrmax,rrmin,zrmin,...
     rzmax,zzmax,rzmin,zzmin] = shapmex(rq,zq,rB,zB,rA,zA,dr2FA,dz2FA,dr,dz) */

# define MEXNAME shapmex

# define RGEOM   pout[ 0]
# define ZGEOM   pout[ 1]
# define AMINO   pout[ 2]
# define EPSIL   pout[ 3]
# define KAPPA   pout[ 4]
# define DELTA   pout[ 5]
# define DELTL   pout[ 6]
# define DELTU   pout[ 7]
# define RRMAX   pout[ 8]
# define ZRMAX   pout[ 9]
# define RRMIN   pout[10]
# define ZRMIN   pout[11]
# define RZMAX   pout[12]
# define ZZMAX   pout[13]
# define RZMIN   pout[14]
# define ZZMIN   pout[15]

# define RQ      prhs[ 0]  /* double, def size=[noq,npq] */
# define ZQ      prhs[ 1]  /* double, size=noq*npq */
# define RB      prhs[ 2]  /* numeric, scalar or empty */
# define ZB      prhs[ 3]  /* numeric, scalar or empty, match size of RB*/
# define RA      prhs[ 4]  /* numeric, scalar */
# define ZA      prhs[ 5]  /* numeric, scalar */
# define DR2FA   prhs[ 6]  /* numeric, scalar */
# define DZ2FA   prhs[ 7]  /* numeric, scalar */
# define DR      prhs[ 8]  /* numeric, scalar */
# define DZ      prhs[ 9]  /* numeric, scalar */
 
 CHECK_NARGIN_EQ(10);

 /* n theta grid, n rho grid, total number of elements */
 size_t noq = mxGetM(RQ), npq = mxGetN(RQ), noo=npq+1,
        nB = mxGetNumberOfElements(RB);

 CHECK_DOUBLE(RQ);
 CHECK_DOUBLE(ZQ);
 CHECK_NUMEL(ZQ,noq*npq);
 CHECK_NUMERIC(RB);
 CHECK_SCALAR_OR_EMPTY(RB);
 CHECK_NUMERIC(ZB);
 CHECK_NUMEL(ZB,nB);
 CHECK_NUMERIC(RA);
 CHECK_SCALAR(RA);
 CHECK_NUMERIC(ZA);
 CHECK_SCALAR(ZA);
 CHECK_NUMERIC(DR2FA);
 CHECK_SCALAR(DR2FA);
 CHECK_NUMERIC(DZ2FA);
 CHECK_SCALAR(DZ2FA);
 CHECK_NUMERIC(DR);
 CHECK_SCALAR(DR);
 CHECK_NUMERIC(DZ);
 CHECK_SCALAR(DZ);

 CHECK_NARGOUT_LE(16);

 mxArray *pout[16];

 double rB, zB; /* boundary points */

#ifdef DEBUG
 mexPrintf("Running shapmex in debug mode\n");
 mexPrintf("noq = %d, npq = %d, noo=%d\n",noq,npq,noo);
#endif

 /* RB,ZB or either both empty or both scalars */
 if (mxIsEmpty(RB)){
  rB=-1; zB=-1; /* placeholders */
#ifdef DEBUG
  mexPrintf("rB,zB empty\n");
#endif
 }
 else
 {
  rB=mxGetScalar(RB);
  zB=mxGetScalar(ZB);
#ifdef DEBUG
  mexPrintf("rB,zB of size 1\n");
#endif
 }


 /* Allocate outputs */
 double *rgeom;
 double *zgeom;
 double *amino;
 double *epsil;
 double *kappa;
 double *delta;
 double *deltl;
 double *deltu;
 double *rrmax; // r at r max
 double *zrmax; // z at r max
 double *rrmin;
 double *zrmin;
 double *rzmax;
 double *zzmax;
 double *rzmin;
 double *zzmin;

 /* Allocate and set pointers to outputs */
 RGEOM = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); rgeom = mxGetPr(RGEOM);
 ZGEOM = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); zgeom = mxGetPr(ZGEOM);
 AMINO = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); amino = mxGetPr(AMINO);
 EPSIL = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); epsil = mxGetPr(EPSIL);
 KAPPA = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); kappa = mxGetPr(KAPPA);
 DELTA = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); delta = mxGetPr(DELTA);
 DELTL = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); deltl = mxGetPr(DELTL);
 DELTU = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); deltu = mxGetPr(DELTU);
 RRMAX = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); rrmax = mxGetPr(RRMAX);
 ZRMAX = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); zrmax = mxGetPr(ZRMAX);
 RRMIN = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); rrmin = mxGetPr(RRMIN);
 ZRMIN = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); zrmin = mxGetPr(ZRMIN);
 RZMAX = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); rzmax = mxGetPr(RZMAX);
 ZZMAX = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); zzmax = mxGetPr(ZZMAX);
 RZMIN = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); rzmin = mxGetPr(RZMIN);
 ZZMIN = mxCreateDoubleMatrix((mwSize)noo,1,mxREAL); zzmin = mxGetPr(ZZMIN);

 /* Special case for values at axis */
 rgeom[0] = mxGetScalar(RA);
 zgeom[0] = mxGetScalar(ZA);
 amino[0] = 0;
 epsil[0] = 0;
 kappa[0] = sqrt(mxGetScalar(DR2FA)/mxGetScalar(DZ2FA));
 delta[0] = 0;
 deltl[0] = 0;
 deltu[0] = 0;
 rrmax[0] = mxGetScalar(RA);
 zrmax[0] = mxGetScalar(ZA);
 rrmin[0] = mxGetScalar(RA);
 zrmin[0] = mxGetScalar(ZA);
 rzmax[0] = mxGetScalar(RA);
 zzmax[0] = mxGetScalar(ZA);
 rzmin[0] = mxGetScalar(RA);
 zzmin[0] = mxGetScalar(ZA);


 /* Call (d)shap function that does all the work - (d)oubles only for mex)*/
 /* These fill indices 2:(npq+1) (1-based) */
 dshap(++rgeom,++zgeom,++amino,++epsil,++kappa,++delta,++deltl,++deltu,
       ++rrmax,++zrmax,++rrmin,++zrmin,++rzmax,++zzmax,++rzmin,++zzmin,
       mxGetPr(RQ),mxGetPr(ZQ),rB,zB,mxGetScalar(DR),mxGetScalar(DZ),npq,noq);

 ASSIGN_PLHS;
}
