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

 /* y = respmex(a,b,c,d) */

# define MEXNAME respmex

# define Y pout[0]

# define A prhs[0] /* double or single, def size=[m,nA] */
# define B prhs[1] /* type match A, def size=[*,nB], size=[nA,nB] */
# define C prhs[2] /* type match A, def size=[*,nC], size=[ m,nC] */
# define D prhs[3] /* type match A, def size=[*,nD], size=[nC,nD] */

 CHECK_NARGIN_EQ(4);

 int nA = mxGetN(A);
 int nB = mxGetN(B);
 int nC = mxGetN(C);
 int nD = mxGetN(D);
 int m = mxGetM(A);

 CHECK_REAL(A);
 CHECK_TYPE_MATCH(B,A);
 CHECK_NROWS(B,nA);
 CHECK_TYPE_MATCH(C,A);
 CHECK_NROWS(C,m);
 CHECK_TYPE_MATCH(D,A);
 CHECK_NROWS(D,nC);

 CHECK_NARGOUT_LE(1);

 mxArray *pout[1];
 
 switch (mxGetClassID(A)) {
  case mxDOUBLE_CLASS: {
   double *a = mxGetPr(A);
   double *b = mxGetPr(B);
   double *c = mxGetPr(C);
   double *d = mxGetPr(D);
   double *y = mxGetPr(Y = mxCreateDoubleMatrix(m, nB+nD, mxREAL));
   dresp(y, a, b, c, d, m, nA, nB, nC, nD);
   break;
  }
  case mxSINGLE_CLASS: {
   float *a = (float *)mxGetData(A);
   float *b = (float *)mxGetData(B);
   float *c = (float *)mxGetData(C);
   float *d = (float *)mxGetData(D);
   float *y = (float *)mxGetData(Y = mxCreateNumericMatrix(m, nB+nD, mxSINGLE_CLASS, mxREAL));
   sresp(y, a, b, c, d, m, nA, nB, nC, nD);
   break;
  }
  /* Other cases impossible due to CHECK_REAL call*/
  }

 ASSIGN_PLHS;
}
