/* nag_1d_quad_guass_1 (d01tac) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 * Mark 7 revised, 2001.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <math.h>
#include <nagd01.h>

#ifdef __cplusplus
extern "C" {
#endif
  static double fun1(double x, Nag_User *comm);
  static double fun2(double x, Nag_User *comm);
  static double fun3(double x, Nag_User *comm);
  static double fun4(double x, Nag_User *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{
  static Integer nstor[3] = {4, 8, 16};

  double a, b;
  Integer i;
  static NagError fail;
  double ans;
  Nag_GaussFormulae gaussformula;
  Nag_Boolean success = Nag_TRUE;
  Nag_User comm;

  fail.print = Nag_TRUE;

  Vprintf("nag_1d_quad_gauss_1 (d01tac) Example Program Results\n");
  Vprintf("\nGauss-Legendre example\n\n");
  for (i=0; i<3; ++i)
    {
      a = 0.0;
      b = 1.0;
      gaussformula = Nag_Legendre;
      /* nag_1d_quad_gauss_1 (d01tac).
       * One-dimensional Gaussian quadrature rule evaluation,
       * thread-safe
       */
      ans = nag_1d_quad_gauss_1(gaussformula, fun1, a, b, nstor[i], &comm,
                                &fail);
      if (fail.code == NE_NOERROR || fail.code == NE_QUAD_GAUSS_NPTS_RULE)
        Vprintf("%ld Points    Answer = %10.5f\n\n", nstor[i], ans);
      else
        {
          Vprintf("%s\n", fail.message);
          success = Nag_FALSE;
        }
    }
  Vprintf("\nGauss-Rational example\n\n");
  for (i=0; i<3; ++i)
    {
      a = 2.0;
      b = 0.0;
      gaussformula = Nag_Rational;
      /* nag_1d_quad_gauss_1 (d01tac), see above. */
      ans = nag_1d_quad_gauss_1(gaussformula, fun2, a, b, nstor[i], &comm,
                                &fail);
      if (fail.code == NE_NOERROR || fail.code == NE_QUAD_GAUSS_NPTS_RULE)
        Vprintf("%ld Points    Answer = %10.5f\n\n", nstor[i], ans);
      else
        {
          Vprintf("%s\n", fail.message);
          success = Nag_FALSE;
        }
    }
  Vprintf("\nGauss-Laguerre example\n\n");
  for (i=0; i<3; ++i)
    {
      a = 2.0;
      b = 1.0;
      gaussformula = Nag_Laguerre;
      /* nag_1d_quad_gauss_1 (d01tac), see above. */
      ans = nag_1d_quad_gauss_1(gaussformula, fun3, a, b, nstor[i], &comm,
                                &fail);
      if (fail.code == NE_NOERROR || fail.code == NE_QUAD_GAUSS_NPTS_RULE)
        Vprintf("%ld Points    Answer = %10.5f\n\n", nstor[i], ans);
      else
        {
          Vprintf("%s\n", fail.message);
          success = Nag_FALSE;
        }
    }
  Vprintf("\nGauss-Hermite  example\n\n");
  for (i=0; i<3; ++i)
    {
      a = -1.0;
      b = 3.0;
      gaussformula = Nag_Hermite;
      /* nag_1d_quad_gauss_1 (d01tac), see above. */
      ans = nag_1d_quad_gauss_1(gaussformula, fun4, a, b, nstor[i], &comm,
                                &fail);
      if (fail.code == NE_NOERROR || fail.code == NE_QUAD_GAUSS_NPTS_RULE)
        Vprintf("%ld Points    Answer = %10.5f\n\n", nstor[i], ans);
      else
        {
          Vprintf("%s\n", fail.message);
          success = Nag_FALSE;
        }
    }
  if (success)
    return EXIT_SUCCESS;
  else
    return EXIT_FAILURE;
}
static double fun1(double x, Nag_User *comm)
{
  return 4.0/(x*x+1.0);
}
static double fun2(double x, Nag_User *comm)
{
  return  1.0/(x*x*log(x));
}
static double fun3(double x, Nag_User *comm)
{
  return  exp(-x)/x;
}
static double fun4(double x, Nag_User *comm)
{
  return  exp(x*(-3.0)*x-x*4.0-1.0);
}