/* Example 5
   =========

   Shows how to implement a Scilab wrapper calling a NAG routine from the 
   Fortran library */

#include "stack-c.h"
#undef Complex
#define Complex NagComplex
#define SciComplex doublecomplex

#include "stdio.h"

// The  header stuff below has been taken from nagmk21_cf_us.h as otherwise
// there are conflicts between doubley defined routines such as dcopy etc.

#ifndef NAG_FTN_INCLUDED
#define NAG_FTN_INCLUDED
#endif

#ifndef RETURN_COMPLEX
#define RETURN_COMPLEX
#endif

#ifdef RETURN_COMPLEX_PARAM
#undef RETURN_COMPLEX
#endif

#ifndef NAG_TYPES
  typedef struct { double re,im; } Complex;

  typedef struct { float re,im; } Complexf;

  typedef int Integer;

  typedef int logical;
#else
  typedef NAG_F77_int logical;
#endif

#define CONST const

#ifdef __cplusplus
#define REFPTR &
#define VOID
#else
#define REFPTR * 
#define VOID void
#endif

#if defined (_WIN32) || defined (_WIN64)
#define NAG_CALL __stdcall
#else
#define NAG_CALL 
#endif

#ifdef __cplusplus
extern "C" {
#endif


#if defined (_WIN32) || defined (_WIN64)
extern void NAG_CALL G05KBF(
  Integer REFPTR igen,
  Integer iseed[] /* 1 dimension */
);
#define \
g05kbf_(igen, iseed)  \
G05KBF(igen, iseed) 
#else
extern void NAG_CALL g05kbf_(
  Integer REFPTR igen,
  Integer iseed[] /* 1 dimension */
);
#endif


#if defined (_WIN32) || defined (_WIN64)
extern void NAG_CALL G05KCF(
  Integer REFPTR igen,
  Integer iseed[] /* 1 dimension */
);
#define \
g05kcf_(igen, iseed)  \
G05KCF(igen, iseed) 
#else
extern void NAG_CALL g05kcf_(
  Integer REFPTR igen,
  Integer iseed[] /* 1 dimension */
);
#endif


#if defined (_WIN32) || defined (_WIN64)
extern void NAG_CALL G05LAF(
  CONST double REFPTR xmu,
  CONST double REFPTR var,
  CONST Integer REFPTR n,
  double x[] /* 1 dimension */,
  CONST Integer REFPTR igen,
  Integer iseed[] /* 1 dimension */,
  Integer REFPTR ifail
);
#define \
g05laf_(xmu, var, n, x, igen, iseed, ifail)  \
G05LAF(xmu, var, n, x, igen, iseed, ifail) 
#else
extern void NAG_CALL g05laf_(
  CONST double REFPTR xmu,
  CONST double REFPTR var,
  CONST Integer REFPTR n,
  double x[] /* 1 dimension */,
  CONST Integer REFPTR igen,
  Integer iseed[] /* 1 dimension */,
  Integer REFPTR ifail
);
#endif

int nag_intext5(char *fname)
{
  // to call this function in scilab use:
  // [x,ifail] = nag_psdrnd_norm_fun(n,xmu,var,repeatable)
  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 n, i, min, max;
  int repeatable;

  Integer igen, ifail, ndim;
  Integer *iseed=0;

  double xmu, var;
  double *x;

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

  GetRhsVar(1, "i", &m1, &n1, &l1); // input n 
  GetRhsVar(2, "d", &m2, &n2, &l2); // input xmu 
  GetRhsVar(3, "d", &m3, &n3, &l3); // input var 
  GetRhsVar(4, MATRIX_OF_BOOLEAN_DATATYPE, &m4, &n4, &l4); // repeatable 

  // there are also other datatypes for integers, real, string etc.
  // MATRIX_OF_INTEGER_DATATYPE on the istk
  // MATRIX_OF_RATIONAL_DATATYPE on the sstk
  // MATRIX OF DOUBLE DATATYPE on the stk
  // MATRIX_OF_COMPLEX_DATATYPE on the zstk
  // MATRIX_OF_BOOLEAN_DATATYPE on the istk
  // STRING_DATATYPE on the cstk 


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

  // Read in n 
  //*============*
  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 output array x 
  // ================================== 
  if ( !(x = malloc(n*sizeof(double))))
    {
      sciprint("Allocation failure\n");
      Error(999); goto END;

    }


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

  // Read in var 
  //==============
  if (m3!=1 || n3!=1)
    {
      sciprint("%s: Dimension should be 1x1 for arg 3\r\n",fname);
      Error(999); goto END;
    }
  else
    {
      var = *stk(l3);
    }

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


  //***********************************
  // End of reading in input variables 
  //***********************************


  if ( !(iseed = malloc(4*sizeof(Integer))))
    {
      sciprint("Allocation failure\n");
      Error(999); goto END;

    }

  // Call initialisation routine for this rng 
  // ======================================== 
  igen = 0; // Basic generator 

  iseed[0] = 23;
  iseed[1] = 4;
  iseed[2] = 10;
  iseed[3] = 1985;

  // Initialise the rng either repeatably or not repeatably
  if (repeatable)
    {
      g05kbf_(&igen, iseed);
    }
  else
    {
      g05kcf_(&igen, iseed);
    }

  // Call random number generator based on Normal PDF
  // * ================================================ 
  ifail = -1;

  g05laf_(&xmu, &var, &ndim, x, &igen, iseed, &ifail);

  // Print NAG error message if routine fails
  //=========================================
  if (ifail != 0)
    sciprint("g05laf failed");

  // Create output variables for Scilab 
  //====================================
  CreateVar(5, "d", &n, &m1, &l5); // x 
  CreateVar(6, "i", &n1, &m1, &l6); // ifail 

  // Assign scilab output variables 
  //================================
  for (i=0;i<n;i++)
    {
      *stk(l5+i) = x[i];
    }

  *istk(l6) = (int) ifail;

  // Assign order for output variables 
  //===================================
  LhsVar(1) = 5;
  LhsVar(2) = 6;

 END:

  // Deallocate memory 
  //===================
  if (x) free(x);
  if (iseed) free(iseed);

  return(0);
}