/* nag_1d_quad_gauss (d01bac) Example Program.
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 * Mark 7 revised, 2001.
 *
 */

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

#ifdef __cplusplus
extern "C" {
#endif
static double NAG_CALL fun1(double x);
static double NAG_CALL fun2(double x);
static double NAG_CALL fun3(double x);
static double NAG_CALL fun4(double x);
#ifdef __cplusplus
}
#endif

int main(int argc, char *argv[])
{
  FILE              *fpout;
  static Integer    nstor[3] = { 4, 8, 16 };
  Integer           exit_status = 0;
  double            a, b;
  Integer           i;
  NagError          fail;
  double            ans;
  Nag_GaussFormulae gaussformula;
  fail.print = Nag_TRUE;

  INIT_FAIL(fail);

  /* Check for command-line IO options */
  fpout = nag_example_file_io(argc, argv, "-results", NULL);
  fprintf(fpout, "nag_1d_quad_gauss (d01bac) Example Program Results\n");
  fprintf(fpout, "\nGauss-Legendre example\n\n");
  for (i = 0; i < 3; ++i)
    {
      a = 0.0;
      b = 1.0;
      gaussformula = Nag_Legendre;
      /* nag_1d_quad_gauss (d01bac).
       * One-dimensional Gaussian quadrature rule evaluation
       */
      ans = nag_1d_quad_gauss(gaussformula, fun1, a, b, nstor[i], &fail);
      if (fail.code == NE_NOERROR || fail.code == NE_QUAD_GAUSS_NPTS_RULE)
        fprintf(fpout, "%ld Points    Answer = %10.5f\n\n", nstor[i], ans);
      else
        {
          fprintf(fpout, "%s\n", fail.message);
          exit_status = 1;
          goto END;
        }
    }
  fprintf(fpout, "\nGauss-Rational example\n\n");
  for (i = 0; i < 3; ++i)
    {
      a = 2.0;
      b = 0.0;
      gaussformula = Nag_Rational;
      /* nag_1d_quad_gauss (d01bac), see above. */
      ans = nag_1d_quad_gauss(gaussformula, fun2, a, b, nstor[i], &fail);
      if (fail.code == NE_NOERROR || fail.code == NE_QUAD_GAUSS_NPTS_RULE)
        fprintf(fpout, "%ld Points    Answer = %10.5f\n\n", nstor[i], ans);
      else
        {
          fprintf(fpout, "%s\n", fail.message);
          exit_status = 1;
          goto END;
        }
    }
  fprintf(fpout, "\nGauss-Laguerre example\n\n");
  for (i = 0; i < 3; ++i)
    {
      a = 2.0;
      b = 1.0;
      gaussformula = Nag_Laguerre;
      /* nag_1d_quad_gauss (d01bac), see above. */
      ans = nag_1d_quad_gauss(gaussformula, fun3, a, b, nstor[i], &fail);
      if (fail.code == NE_NOERROR || fail.code == NE_QUAD_GAUSS_NPTS_RULE)
        fprintf(fpout, "%ld Points    Answer = %10.5f\n\n", nstor[i], ans);
      else
        {
          fprintf(fpout, "%s\n", fail.message);
          exit_status = 1;
          goto END;
        }
    }
  fprintf(fpout, "\nGauss-Hermite  example\n\n");
  for (i = 0; i < 3; ++i)
    {
      a = -1.0;
      b = 3.0;
      gaussformula = Nag_Hermite;
      /* nag_1d_quad_gauss (d01bac), see above. */
      ans = nag_1d_quad_gauss(gaussformula, fun4, a, b, nstor[i], &fail);
      if (fail.code == NE_NOERROR || fail.code == NE_QUAD_GAUSS_NPTS_RULE)
        fprintf(fpout, "%ld Points    Answer = %10.5f\n\n", nstor[i], ans);
      else
        {
          fprintf(fpout, "%s\n", fail.message);
          exit_status = 1;
          goto END;
        }
    }

 END:
  if (fpout != stdout) fclose(fpout);
  return exit_status;
}


static double NAG_CALL fun1(double x)
{
  return 4.0/(x*x+1.0);
}

static double NAG_CALL fun2(double x)
{
  return 1.0/(x*x*log(x));
}

static double NAG_CALL fun3(double x)
{
  return exp(-x)/x;
}

static double NAG_CALL fun4(double x)
{
  return exp(x*(-3.0)*x-x*4.0-1.0);
}