/* [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved. */
# include "genlib.h"

void mexFunction(int nlhs,       mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{

# define MEXNAME greenemmex

 /* [y,...] = greenemmex(r1,z1,r2,z2,mode); */

# define R1     prhs[0] /* double, def size=n1 */
# define Z1     prhs[1] /* double, size=n1 */
# define R2     prhs[2] /* double, def size=n2 */
# define Z2     prhs[3] /* double, size=n1 */
# define MODE   prhs[4] /* int32  */

 CHECK_NARGIN_EQ(5);

 int n1 = mxGetNumberOfElements(R1), n2 = mxGetNumberOfElements(R2),
     nm = mxGetNumberOfElements(MODE);

 CHECK_DOUBLE(R1);
 CHECK_DOUBLE(Z1);
 CHECK_NUMEL(Z1,n1);
 CHECK_DOUBLE(R2);
 CHECK_DOUBLE(Z2);
 CHECK_NUMEL(Z2,n2);
 CHECK_INT32(MODE);

 CHECK_NARGOUT_LE(nm);

 bool der1 = 0, der2 = 0;
 int i,j,k;
 int    *mode = (int *)mxGetData(MODE);
 double **y = mxCalloc(nlhs, sizeof(double *)),
       *pr1 = mxGetPr(R1), *pr2 = mxGetPr(R2), 
       *pz1 = mxGetPr(Z1), *pz2 = mxGetPr(Z2);

 for (k=0; k<nlhs; k++) {
  y[k] = mxGetPr(plhs[k] = mxCreateDoubleMatrix(n1,n2,mxREAL));
  der1 |= mode[k] >  0;
  der2 |= mode[k] > 10;
 }
 
 for (j=0; j<n2; j++) {
  for (i=0; i<n1; i++) {
   
   greenem(mode, nm, der1, der2, pr1[i], pz1[i], pr2[j], pz2[j], y);

   for (k=0; k<nm; k++) {
    y[k]++;
   } /* for k */
  } /* for i */
 } /* for j */
}
