/* nag_quad_1d_gauss_wrec (d01tdc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.2, 2017.
 */

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

int main(void)
{
  /* Scalars */
  Integer      exit_status = 0;
  Integer      n, i;
  double       muzero, alpha, beta, a1, b1, ab1, ab2, d, ri, fa1, fb1, fab2;
  /* Arrays */
  char         nag_enum_arg[40];
  double       *a = 0, *b = 0, *c = 0, *abscis = 0, *weight = 0;
  const char   *str_type;
  /* Nag Types */
  Nag_QuadType quadtype;
  NagError     fail, fail1, fail2;

  INIT_FAIL(fail);

  printf("nag_quad_1d_gauss_wrec (d01tdc) Example Program Results\n");
  /* Skip heading in data file */
  scanf("%*[^\n] ");
  /* Input number of abscissae required, n */
  scanf("%" NAG_IFMT "%*[^\n] ", &n);
  /* Input quadrature rule to simulate, Nag_QuadType */
  scanf("%39s%*[^\n] ", nag_enum_arg);
  quadtype = (Nag_QuadType) nag_enum_name_to_value(nag_enum_arg);

  /* Allocate coefficient, weight and abscissae arrays */
  if (!(a = NAG_ALLOC(n, double)) ||
      !(b = NAG_ALLOC(n, double)) ||
      !(c = NAG_ALLOC(n, double)) ||
      !(weight = NAG_ALLOC(n, double)) ||
      !(abscis = NAG_ALLOC(n, double)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Convert QuadType to string using
   * nag_enum_value_to_name (x04nbc).
   * Converts NAG enum member value to its name
   */
  str_type = nag_enum_value_to_name(quadtype);
  printf("\nQuadrature type = %s\n\n",str_type);

  /* Generate recursion coefficients a, b, c from quadtype */
  switch (quadtype) {
  case Nag_Quad_Gauss_Legendre:
    {
      muzero = 2.0;
      for (i = 0; i < n; ++i) {
        ri = (double) (i);
        a[i] = (2.0*ri + 1.0)/(ri + 1.0);
        b[i] = 0.0;
        c[i] = ri/(ri + 1.0);
      }
    }
    break;
  case Nag_Quad_Gauss_Chebyshev_first:
    {
      muzero = X01AAC;
      for (i = 0; i < n; ++i) {
        a[i] = 2.0;
        b[i] = 0.0;
        c[i] = 1.0;
      }
      a[0] = 1.0;
    }
    break;
  case Nag_Quad_Gauss_Chebyshev_second:
    {
      muzero = 0.5*X01AAC;
      for (i = 0; i < n; ++i) {
        a[i] = 2.0;
        b[i] = 0.0;
        c[i] = 1.0;
      }
    }
    break;
  case Nag_Quad_Gauss_Jacobi:
    {
      /* Input quadrature paramaters alpha and beta */
      scanf("%lf%lf\n", &alpha, &beta);
      printf("Using parameters alpha = %10.5f, beta = %10.5f\n\n", alpha, beta);
      a1   = alpha + 1.0;
      b1   = beta  + 1.0;
      ab1  = alpha + beta + 1.0;
      ab2  = a1    + b1;
      INIT_FAIL(fail1);
      INIT_FAIL(fail2);
      /* nag_gamma (s14aac).
       * Gamma function Gamma(x)
       */
      fa1  = nag_gamma(a1, &fail);
      fb1  = nag_gamma(b1, &fail1);
      fab2 = nag_gamma(ab2, &fail2);
      if (fail.code != NE_NOERROR || fail1.code != NE_NOERROR ||
          fail2.code != NE_NOERROR) {
        printf("Error from nag_gamma (s14aac).\n%s\n", fail.message);
        exit_status = 1;
        goto END;
      }
      muzero = pow(2.0,ab1)*fa1*fb1/fab2;
      a[0] = 0.5*ab2;
      b[0] = 0.5*(alpha-beta);
      c[0] = 0.0;
      for (i = 1; i < n; ++i) {
        ri = (double) i;
        d    = (2.0*ri + 2.0)*(ri + ab1);
        a[i] = (2.0*ri + ab1)*(2.0*ri + ab2)/d;
        d    = (2.0*ri + alpha + beta)*d;
        b[i] = (2.0*ri + ab1)*(alpha*alpha - beta*beta)/d;
        c[i] = 2.0*(ri + alpha)*(ri + beta)*(2.0*ri + ab2)/d;
      }
    }
    break;
  case Nag_Quad_Gauss_Laguerre:
    {
      /* Input quadrature paramaters alpha */
      scanf("%lf\n", &alpha);
      printf("Using parameter alpha = %10.5f\n\n",alpha);
      a1   = alpha + 1.0;
      /* nag_gamma (s14aac).
       * gamma(x)
       */
      muzero  = nag_gamma(a1, &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_gamma (s14aac).\n%s\n", fail.message);
        exit_status = 1;
        goto END;
      }
      for (i = 0; i < n; ++i) {
        ri = (double) i;
        a[i] = -1.0/(ri + 1.0);
        b[i] = (2.0*ri + a1)/(ri + 1.0);
        c[i] = (ri + alpha)/(ri + 1.0);
      }
    }
    break;
  case Nag_Quad_Gauss_Hermite:
    {
      muzero = sqrt(X01AAC);
      for (i = 0; i < n; ++i) {
        a[i] = 2.0;
        b[i] = 0.0;
        c[i] = (double) 2*i;
      }
    }
    break;
  default:
    {
      exit_status = 1;

      printf("The quadrature type %s is not handled here\n", str_type);
      goto END;
    }
  }

  /* nag_quad_1d_gauss_wrec (d01tdc).
   * Compute weights and abscissae for a Gaussian quadrature rule
   * governed by a three-term recurrence relation.
   */
  nag_quad_1d_gauss_wrec(n, a, b, c, muzero, weight, abscis, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_quad_1d_gauss_wrec (d01tdc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  printf(" %3" NAG_IFMT " points\n\n    Abscissae        Weights\n\n", n);
  for (i = 0; i < n; i++) {
    printf("%15.6e", abscis[i]);
    printf("%15.6e\n", weight[i]);
  }
  printf("\n");

END:
  NAG_FREE(a);
  NAG_FREE(b);
  NAG_FREE(c);
  NAG_FREE(abscis);
  NAG_FREE(weight);

  return exit_status;
}