/* Example 4b
   =========

   Shows how to implement a Scilab wrapper calling a NAG routine using a User
   defined function written in Scilab.  This requires an special C interface
   for this second function */

#include "stack-c.h"
#undef Complex
#define Complex NagComplex
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd01.h>
#define SciComplex doublecomplex

#include "stdio.h"

static double fSci(Integer n, double z[], Nag_User *comm);

int nag_intext4b(char *fname)
{
  /* to call this function in scilab use:
   [finval,acc,minpts,fail] = nag_multid_quad_funb(ndim,a,b,minpts,maxpts,eps)*/
  int m1,n1,l1;
  int m2,n2,l2;
  int m3,n3,l3;
  int m4,n4,l4;
  int m5,n5,l5;
  int m6,n6,l6;
  int l7;
  int l8;
  int l9;
  int l10;

  int n, i, min, max;

  Integer ndim, minpts, maxpts;
  double eps, finval, acc;
  double *a=0, *b=0;
  Nag_User comm;
  NagError ifail;

  /* define minimum and maximum left and right hand arguments*/
  int minlhs=1, minrhs=6, maxlhs=4, maxrhs=6;
  Nbvars = 0;
  CheckRhs(minrhs, maxrhs);
  CheckLhs(minlhs,maxlhs);

/* Initialise ifail */
  INIT_FAIL(ifail);

  GetRhsVar(1, "i", &m1, &n1, &l1); /* input ndim */
  GetRhsVar(2, "d", &m2, &n2, &l2); /* input a */
  GetRhsVar(3, "d", &m3, &n3, &l3); /* input b */
  GetRhsVar(4, "i", &m4, &n4, &l4); /* input minpts */
  GetRhsVar(5, "i", &m5, &n5, &l5); /* input maxpts */
  GetRhsVar(6, "d", &m6, &n6, &l6); /* input eps */


/********************************************************************/
/* Read in input variables from Scilab and convert to NAG variables */
/********************************************************************/

/* Read in ndim */
/**============**/
  if (m1!=1 || n1!=1)
    {
      sciprint("%s: Dimension should be 1x1 character for arg 1\r\n",fname);
      Error(999); goto END;
    }
  else
    {
      n = *istk(l1);
      ndim = (Integer) n;
    }

  /* Allocate memory for a(ndim) and b(ndim)*/
  /*========================================*/
  if ( !(a = NAG_ALLOC(n,double) ) ||
       !(b = NAG_ALLOC(n,double)))
    {
      sciprint("Allocation failure\n");
      Error(999); goto END;
    }


  /* Read in a */
  /*===========*/
  if (m2!= n || n2!=1)
    {
      sciprint("%s: Dimension should be ndim by 1 for arg 2\r\n",fname);
      sciprint("n = %i", n);
      Error(999); goto END;
    }
  else
    {
      for (i=0; i<n; i++)
        {
          a[i] = *stk(l2+i);
        }
    }

  /* Read in b */
  /*==============*/
  if (m3!=n || n3!=1)
    {
      sciprint("%s: Dimension should be ndim by 1 for arg 3\r\n",fname);
      Error(999); goto END;
    }
  else
    {
      for(i=0; i<n; i++)
        {
          b[i] = *stk(l3+i);
        }
    }

  /* Read in minpts */
  /*============*/
  if (m4!=1 || n4!=1)
    {
      sciprint("%s: Dimension should be 1x1 for arg 4\r\n",fname);
      Error(999); goto END;
    }
  else
    {
      min = *istk(l4);
      minpts = (Integer)min;
    }

  /* Read in maxpts */
  /*============*/
  if (m5!=1 || n5!=1)
    {
      sciprint("%s: Dimension should be 1x1 for arg 5\r\n",fname);
      Error(999); goto END;
    }
  else
    {
      max = *istk(l5);
      maxpts = (Integer) max;
    }

  /* Read in eps */
  /*==============*/
  if (m6!=1 || n6!=1)
    {
      sciprint("%s: Dimension should be 1x1 for arg 6\r\n",fname);
      Error(999); goto END;
    }
  else
    {
      eps = *stk(l6);
    }
  /*************************************/
  /* End of reading in input variables */
  /*************************************/

  /* nag_multid_quad_adapt_1 (d01wcc).
   * Multi-dimensional adaptive quadrature, thread-safe
   *====================================================*/
  nag_multid_quad_adapt_1(ndim, fSci, a, b, &minpts, maxpts, eps, &finval, &acc,
                          &comm, &ifail);

  /* Print NAG error message if routine fails*/
  /*=========================================*/
  if (ifail.code != NE_NOERROR)
    sciprint("Ifail message is %s", ifail.message);

  /* Create output variables for Scilab */
  /*====================================*/
  CreateVar(7, "d", &n1, &m1, &l7); /* finval */
  CreateVar(8, "d", &n1, &m1, &l8); /* acc */
  CreateVar(9, "i", &n1, &m1, &l9); /* minpts */
  CreateVar(10, "i", &n1, &m1, &l10); /* ifail */

  /* Assign scilab output variables */
  /*================================*/
  *stk(l7) = finval;
  *stk(l8) = acc;
  *istk(l9) = (int) minpts;
  *istk(l10) = (int) ifail.code;

  /* Assign order for output variables */
  /*===================================*/
  LhsVar(1) = 7;
  LhsVar(2) = 8;
  LhsVar(3) = 9;
  LhsVar(4) = 10;


 END:

  /* Deallocate memory */
  /*===================*/
  if (a) NAG_FREE(a);
  if (b) NAG_FREE(b);

  return(0);
}

static double fSci(Integer n, double z[], Nag_User *comm)
{
  int mlhs, mrhs, ibegin;
  int m_u, n_u, l_u;
  int m_z, n_z, l_z;
  int i;
  double retval;
  char name[] = "fUser";

  // Put input variables for scilab function onto stack in the right order  
  // On input u (stored in space 7) is a dummy input to be overwritten with 
  // the output from fUser 
  // The inputs for fUser must be in the correct order in the last places of 
  // the stack.  They are then scrapped and the outputs are created starting 
  // at the same place as the first input variable.
  m_u = 1;
  n_u = 1;
  m_z = (int) n;
  n_z = 1;
  
  CreateVar(7, "d", &m_u, &n_u, &l_u); // first input 
  CreateVar(8, "d", &m_z, &n_z, &l_z); // second input 

  // Write the input vector z onto the stack in the second input 
  for (i=0; i<m_z; i++)
    {
      *stk(l_z+i) = z[i];
    }


  ibegin = 7; // the first input variable on the stack
  mlhs = 1;  // number of left hand side outputs to write onto the stack
  mrhs = 2; // number of right hand side inputs to read off the stack
  SciString(&ibegin,name,&mlhs,&mrhs); // Call scilab function 'name'

  // output from scilab function is now on the stack so return the result 
  retval = *stk(l_u);

  return retval;
}