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

 /* [xa,ya,fa,dx2fa,dy2fa,dxyfa,xs,ys,fs,dx2fs,dy2fs,dxyfs,stat] = asxymex(fxy,x,y,dc,dx,dy,idx,idy,lxy[,dimw]) */

# define MEXNAME asxymex

# define XA    pout[ 0]
# define YA    pout[ 1]
# define FA    pout[ 2]
# define DX2FA pout[ 3]
# define DY2FA pout[ 4]
# define DXYFA pout[ 5]
# define XS    pout[ 6]
# define YS    pout[ 7]
# define FS    pout[ 8]
# define DX2FS pout[ 9]
# define DY2FS pout[10]
# define DXYFS pout[11]
# define STAT  pout[12]

# define FXY   prhs[ 0]  /* double, size=nx*ny        */
# define X     prhs[ 1]  /* double, def size=nx       */
# define Y     prhs[ 2]  /* double, def size=ny       */
# define DC    prhs[ 3]  /* numeric, scalar           */
# define DX    prhs[ 4]  /* numeric, scalar           */
# define DY    prhs[ 5]  /* numeric, scalar           */
# define IDX   prhs[ 6]  /* numeric, scalar           */
# define IDY   prhs[ 7]  /* numeric, scalar           */
# define LXY   prhs[ 8]  /* logical, size=nx*ny       */
# define DIMW  prhs[ 9]  /* numeric, scalar, optional */

 CHECK_NARGIN_GE( 9);
 CHECK_NARGIN_LE(10);
 
 int nx = mxGetNumberOfElements(X), ny = mxGetNumberOfElements(Y),
     dima, na, dims, ns;

 CHECK_DOUBLE(FXY);
 CHECK_NUMEL(FXY,nx*ny);
 CHECK_DOUBLE(X);
 CHECK_DOUBLE(Y);
 CHECK_NUMERIC(DC);
 CHECK_SCALAR(DC);
 CHECK_NUMERIC(DX);
 CHECK_SCALAR(DX);
 CHECK_NUMERIC(DY);
 CHECK_SCALAR(DY);
 CHECK_NUMERIC(IDX);
 CHECK_SCALAR(IDX);
 CHECK_NUMERIC(IDY);
 CHECK_SCALAR(IDY);
 CHECK_LOGICAL(LXY);
 CHECK_NUMEL(LXY,nx*ny);
 if (nrhs == 10) { CHECK_NUMERIC(DIMW);
                   CHECK_SCALAR(DIMW);  }

 CHECK_NARGOUT_LE(13);

 mxArray *pout[13];

 if (nrhs == 10) dima = dims = mxGetScalar(DIMW);
 else            dima = dims = nx*ny;
 
 double *xa    = (double *)mxCalloc(  dima,sizeof(double));
 double *ya    = (double *)mxCalloc(  dima,sizeof(double));
 double *fa    = (double *)mxCalloc(  dima,sizeof(double));
 double *dx2fa = (double *)mxCalloc(  dima,sizeof(double));
 double *dy2fa = (double *)mxCalloc(  dima,sizeof(double));
 double *dxyfa = (double *)mxCalloc(  dima,sizeof(double));
 double *xs    = (double *)mxCalloc(  dims,sizeof(double));
 double *ys    = (double *)mxCalloc(  dims,sizeof(double));
 double *fs    = (double *)mxCalloc(  dims,sizeof(double));
 double *dx2fs = (double *)mxCalloc(  dims,sizeof(double));
 double *dy2fs = (double *)mxCalloc(  dims,sizeof(double));
 double *dxyfs = (double *)mxCalloc(  dims,sizeof(double));
 int    *work  = (int    *)mxCalloc(3*dims,sizeof(int   ));
 STAT=mxCreateLogicalScalar((mxLogical) dasxy(xa,ya,fa,dx2fa,dy2fa,dxyfa,dima,&na,
                                              xs,ys,fs,dx2fs,dy2fs,dxyfs,dims,&ns,
                                              mxGetPr(FXY),FLTC(0.0),mxGetPr(X),mxGetPr(Y),mxGetScalar(DC),nx,ny,
                                              mxGetScalar(DX),mxGetScalar(DY),mxGetScalar(IDX),mxGetScalar(IDY),
                                              (bool *)mxGetData(LXY),work));
 memcpy((void *)mxGetPr(XA    = mxCreateDoubleMatrix(na,1,mxREAL)),(void *)xa   ,na*sizeof(double));
 memcpy((void *)mxGetPr(YA    = mxCreateDoubleMatrix(na,1,mxREAL)),(void *)ya   ,na*sizeof(double));
 memcpy((void *)mxGetPr(FA    = mxCreateDoubleMatrix(na,1,mxREAL)),(void *)fa   ,na*sizeof(double));
 memcpy((void *)mxGetPr(DX2FA = mxCreateDoubleMatrix(na,1,mxREAL)),(void *)dx2fa,na*sizeof(double));
 memcpy((void *)mxGetPr(DY2FA = mxCreateDoubleMatrix(na,1,mxREAL)),(void *)dy2fa,na*sizeof(double));
 memcpy((void *)mxGetPr(DXYFA = mxCreateDoubleMatrix(na,1,mxREAL)),(void *)dxyfa,na*sizeof(double));
 memcpy((void *)mxGetPr(XS    = mxCreateDoubleMatrix(ns,1,mxREAL)),(void *)xs   ,ns*sizeof(double));
 memcpy((void *)mxGetPr(YS    = mxCreateDoubleMatrix(ns,1,mxREAL)),(void *)ys   ,ns*sizeof(double));
 memcpy((void *)mxGetPr(FS    = mxCreateDoubleMatrix(ns,1,mxREAL)),(void *)fs   ,ns*sizeof(double));
 memcpy((void *)mxGetPr(DX2FS = mxCreateDoubleMatrix(ns,1,mxREAL)),(void *)dx2fs,ns*sizeof(double));
 memcpy((void *)mxGetPr(DY2FS = mxCreateDoubleMatrix(ns,1,mxREAL)),(void *)dy2fs,ns*sizeof(double));
 memcpy((void *)mxGetPr(DXYFS = mxCreateDoubleMatrix(ns,1,mxREAL)),(void *)dxyfs,ns*sizeof(double));
 mxFree(   xa);
 mxFree(   ya);
 mxFree(   fa);
 mxFree(dx2fa);
 mxFree(dy2fa);
 mxFree(dxyfa);
 mxFree(   xs);
 mxFree(   ys);
 mxFree(   fs);
 mxFree(dx2fs);
 mxFree(dy2fs);
 mxFree(dxyfs);
 mxFree( work);

 ASSIGN_PLHS;
}
