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

 /* Fx = gszrmex(Fb,Jy,cx,cq,cr,cs,ci,co,dz) */

# define MEXNAME gszrmex

# define FX pout[0]

# define FB prhs[0] /* double, size=2*(nz-nr-2) */
# define JY prhs[1] /* double, def size=[nz-2,nr-2] */
# define CX prhs[2] /* double, size=(nr-2) */
# define CQ prhs[3] /* double, size=(nr-2)*(nz-2) */
# define CR prhs[4] /* double, size=(nr-2)*(nz-2) */
# define CS prhs[5] /* double, size=(nr-2)*(nz-2) */
# define CI prhs[6] /* numeric, scalar */
# define CO prhs[7] /* numeric, scalar */
# define DZ prhs[8] /* numeric, scalar */

  CHECK_NARGIN_EQ(9);
  
  int nz = (int)mxGetM(JY)+2;
  int nr = (int)mxGetN(JY)+2;
  
  CHECK_DOUBLE(FB);
  CHECK_NUMEL(FB,2*(nz+nr-2));
  CHECK_DOUBLE(JY);
  CHECK_DOUBLE(CX);
  CHECK_NUMEL(CX,(nr-2));
  CHECK_DOUBLE(CQ);
  CHECK_NUMEL(CQ,(nr-2)*(nz-2));
  CHECK_DOUBLE(CR);
  CHECK_NUMEL(CR,(nr-2)*(nz-2));
  CHECK_DOUBLE(CS);
  CHECK_NUMEL(CS,(nr-2)*(nz-2));
  CHECK_NUMERIC(CI);
  CHECK_SCALAR(CI);
  CHECK_NUMERIC(CO);
  CHECK_SCALAR(CO);
  CHECK_NUMERIC(DZ);
  CHECK_SCALAR(DZ);

  CHECK_NARGOUT_LE(1);

 mxArray *pout[1];

  double *fx = mxGetPr(FX = mxCreateDoubleMatrix(nz,nr,mxREAL));
  double *fb = mxGetPr    (FB);
  double *jy = mxGetPr    (JY);
  double *cx = mxGetPr    (CX);
  double *cq = mxGetPr    (CQ);
  double *cr = mxGetPr    (CR);
  double *cs = mxGetPr    (CS);
  double  ci = mxGetScalar(CI);
  double  co = mxGetScalar(CO);
  double  dz = mxGetScalar(DZ);
  /* allocate nr-2 x nz-2 for parallelised gszr */
  double *p = (double *)mxCalloc(nr-2,sizeof(double));
  dgszr(fx,fb,jy,cx,cq,cr,cs,ci,co,p,nz,nr,dz);

 ASSIGN_PLHS;
}
