/* The easiest way to pull in all necessary definitions is to include oct.h
   Before we can do it, we need to overwrite Complex type duplicated 
   in Octave and NAG header files */
#undef Complex
#define Complex OctComplex
#include <octave/oct.h>
/* For NAG definitions we need to include the appropriate NAG header files.
   Again, Complex type has to be overwritten first*/
#undef Complex
#define Complex NagComplex
// MatrixType is another duplicated definition
#define MatrixType NagMatrixType
#include <nag.h>
#include <nags.h>

/* The third parameter (nargout) is not used, so it is
   omitted from the list of arguments to DEFUN_DLD in order to avoid
   the warning from gcc about an unused function parameter. */
DEFUN_DLD (nag_cmplx_erfc, args, ,
  "Calls nag_complex_erfc (s15ddc), which computes values \n\
of the function w(z)=e^(-z^2)erfc-iz, for complex z.\n")
{
  // Variable to store function output values
  octave_value_list retval;
  int nargin = args.length ();
  if (nargin != 1)
    {
      error ("Insufficient input arguments.");
    }
  else
    {
      // Retrieve input arguments from args
      OctComplex z = args(0).complex_value();
      if (! error_state)
        {
          // Declare local variables
          NagError fail;
          INIT_FAIL(fail);
          OctComplex w;
          NagComplex z_nag, w_nag;
          /* NAG and Octave Complex types differ, so we need to
             copy the input OctComplex variable z into NagComplex z_nag */
          z_nag.re = z.real();
          z_nag.im = z.imag();
          
          // Call NAG routine like in C
          /* nag_complex_erfc (s15ddc).
           * Scaled complex complement of error function,
           * e^(-z^2)erfc(-iz)
           */
          w_nag = nag_complex_erfc(z_nag,&fail);
          
          if (fail.code!=NE_NOERROR)
            {
              error ("nag_complex_erfc (s15ddc).\n%s\n", fail.message);
              goto END;
            }
          /* The result is store in NAG Complex type, so 
             we need to copy it to Octave Complex type. */
          w.real() = w_nag.re;
          w.imag() = w_nag.im;
          // Assign output arguments to retval
          retval(0) = w;
          retval(1) = fail.code;
        }
    }

 END:
  return retval;
}