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

 /*  [a,dan,r,z,ir] = rtcimex(ai,ero,ezo,Fx,Fq,cr,cz,Opy,Fo,Opo,dap,ro,aro,zo,azo) */

# define MEXNAME rtcimex

# define A   pout[ 0]
# define DAN pout[ 1]
# define R   pout[ 2]
# define Z   pout[ 3]
# define IR  pout[ 4]

# define AI  prhs[ 0]  /* double, def size=[no,np] */
# define ERO prhs[ 1]  /* double, def size=mo */
# define EZO prhs[ 2]  /* double, size=mo */
# define F   prhs[ 3]  /* double, def size=[nz,nr] */
# define FP  prhs[ 4]  /* double, size=np */
# define CR  prhs[ 5]  /* double, size=no */
# define CZ  prhs[ 6]  /* double, size=no */
# define OPY prhs[ 7]  /* int8, size=(nz-2)*(nr-2) */
# define FO  prhs[ 8]  /* double, size=mo */
# define OPO prhs[ 9]  /* int8, size=1 */
# define DAP prhs[10]  /* numeric, scalar */
# define RO  prhs[11]  /* double, size=mo, optional */
# define ARO prhs[12]  /* double, size=no, optional */
# define ZO  prhs[13]  /* double, size=mo, optional */
# define AZO prhs[14]  /* double, size=no, optional */

 /* Replacement for CHECK_NARGIN_* */
 if ((nrhs != 11) && (nrhs != 13) && (nrhs != 15)) mexErrMsgIdAndTxt("rtcimex:nargin","rtcimex should have 11 or 13 or 15 arguments");

 int nz = mxGetM(F), nr=mxGetN(F), no = mxGetM(AI), np = mxGetN(AI), mo = mxGetNumberOfElements(ERO);

 CHECK_DOUBLE(AI);
 CHECK_DOUBLE(ERO);
 CHECK_DOUBLE(EZO);
 CHECK_NUMEL(EZO,mo);
 CHECK_DOUBLE(F);
 CHECK_DOUBLE(FP);
 CHECK_NUMEL(FP,np);
 CHECK_DOUBLE(CR);
 CHECK_NUMEL(CR,no);
 CHECK_DOUBLE(CZ);
 CHECK_NUMEL(CZ,no);
 CHECK_INT8(OPY);
 CHECK_NUMEL(OPY,(nz-2)*(nr-2));
 CHECK_DOUBLE(FO);
 CHECK_NUMEL(FO,mo);
 CHECK_INT8(OPO);
 CHECK_NUMEL(OPO,1);
 CHECK_NUMERIC(DAP);
 CHECK_SCALAR(DAP);
 if (nrhs>11)
 {
  CHECK_DOUBLE(RO);
  CHECK_NUMEL(RO,mo);
  CHECK_DOUBLE(ARO);
  CHECK_NUMEL(ARO,no);
 }
 if (nrhs>13)
 {
  CHECK_DOUBLE(ZO);
  CHECK_NUMEL(ZO,mo);
  CHECK_DOUBLE(AZO);
  CHECK_NUMEL(AZO,no);
 }
 
 if (nrhs == 11) {CHECK_NARGOUT_LE(2);}
 else if (nrhs == 13) {CHECK_NARGOUT_LE(3);}
 else {CHECK_NARGOUT_LE(5);}

 mxArray *pout[5];
 
 int8_t *Opo = (int8_t *) mxGetData(OPO);
 A   = mxCreateDoubleMatrix(no,np,mxREAL);
 R   = mxCreateDoubleMatrix(no,np,mxREAL);
 IR  = mxCreateDoubleMatrix(no,np,mxREAL);
 DAN = mxCreateDoubleMatrix( 1, 1,mxREAL);
 drtci(mxGetPr(A ),mxGetPr(DAN),
       mxGetPr(AI),mxGetPr(ERO),mxGetPr(EZO),
       mxGetPr(F ),mxGetPr(FP ),
       mxGetPr(CR),mxGetPr(CZ ),(int8_t *) mxGetData(OPY),
       mxGetPr(FO),        *Opo,mxGetScalar(DAP),
       nr,nz,np,no,mo);
 if (nlhs > 2) {
  R = mxCreateDoubleMatrix(no,np,mxREAL);
  drtci1(mxGetPr(R),mxGetPr(A),mxGetPr(ARO),mxGetPr(RO),np,no,mo);
 }
 if (nlhs > 3) {
  Z = mxCreateDoubleMatrix(no,np,mxREAL);
  drtci1(mxGetPr(Z),mxGetPr(A),mxGetPr(AZO),mxGetPr(ZO),np,no,mo);
 }
 if (nlhs > 4) {
  IR = mxCreateDoubleMatrix(no,np,mxREAL);
  vdInv(no*np,mxGetPr(R),mxGetPr(IR));
 }

 ASSIGN_PLHS;
}
