/* Example 2
   =========

   Shows how to implement a simple NAG routine. Using complex types
   and variable left hand side arguments is demonstrated. */

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

int nag_intext2(char *fname)
{

  // y = nag_elliptic_integral_f_eqn(z,akp,a,b) 
  int m1,n1,l1;
  int m2,n2,l2;
  int m3,n3,l3;
  int m4,n4,l4;
  int m5,n5,l5;
  int l6;

  double akp, a, b;
  SciComplex z, y;
  NagComplex z_nag, y_nag;
  NagError ifail;
  // 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);

  // Initialise ifail 
  INIT_FAIL(ifail);

  // Get the right hand side variables 1 to 4
  GetRhsVar(1, "z", &m1, &n1, &l1);
  GetRhsVar(2, "d", &m2, &n2, &l2);
  GetRhsVar(3, "d", &m3, &n3, &l3);
  GetRhsVar(4, "d", &m4, &n4, &l4);

//   d stands for double, r for float (real), i for integer and 
//     c for character string and z for complex 

  if (m1!=1 || n1!=1)
    {
      sciprint("%s: Dimension should be 1x1 for arg 1\r\n",fname);
      Error(999); return(0);
    }

  if (m2!=1 || n2!=1)
    {
      sciprint("%s: Dimension should be 1x1 for arg 2\r\n",fname);
      Error(999); return(0);
    }

  if (m3!=1 || n3!=1)
    {
      sciprint("%s: Dimension should be 1x1 for arg 3\r\n",fname);
      Error(999); return(0);
    }

  if (m4!=1 || n4!=1)
    {
      sciprint("%s: Dimension should be 1x1 for arg 4\r\n",fname);
      Error(999); return(0);
    }

  // Create the 5th and 6th variable for the left hand side output variables
  CreateVar(5, "i", &m1, &n1, &l5);
  CreateVar(6, "z", &m1, &n1, &l6);

  // Retrieve the input data from the various stacks. zstk() is for complex, stk()
  // is for double, sstk() is for reals, cstk() is for characters and istk() is
  // for integers.
  z = *zstk(l1);
  akp = *stk(l2);
  a = *stk(l3);
  b = *stk(l4);

  // Convert Scilab to NAG complex types 
  z_nag.re = z.r;
  z_nag.im = z.i;

  // Run the NAG function s21dac 
  y_nag = nag_elliptic_integral_f(z_nag,akp,a,b,&ifail);

  // Convert NAG to Scilab complex types
  y.r = y_nag.re;
  y.i = y_nag.im;

  // Put the outputs of the NAG function into the created variables 5 and 6
  *zstk(l6) = y;
  *istk(l5) = ifail.code;

  // Point the ouput of the Scilab function to the variables 5 and 6. The order matters 
  // here, as we want 6 to be the first variable out, and 5 to be the second optional 
  // output variable.
  LhsVar(1) = 6;
  LhsVar(2) = 5;
  return(0);
}