/* nag_ode_bvp_ps_lin_quad_weights (d02uyc) Example Program.
 *
 * Copyright 2011, Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd02.h>
#include <nagx02.h>

#ifdef __cplusplus
extern "C" {
#endif
  static double NAG_CALL exact(double x);
#ifdef __cplusplus
}
#endif

int main(void)
{
  /*  Scalars */
  Integer     exit_status = 0;
  Integer     i, n;
  double      a = -1.0, b = 3.0;
  double      integ, scale, uerr;
  double      teneps = 10.0 * nag_machine_precision;
  /*  Arrays */
  double      *f = 0, *w = 0, *x = 0;
  /* NAG types */
  Nag_Boolean reqerr = Nag_FALSE, reqwgt = Nag_FALSE;
  NagError    fail;

  INIT_FAIL(fail);

  printf("nag_ode_bvp_ps_lin_quad_weights (d02uyc) "
         "Example Program Results\n\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  scanf("%"NAG_IFMT "%*[^\n] ", &n);
  if (
    !(f = NAG_ALLOC((n + 1), double)) ||
    !(w = NAG_ALLOC((n + 1), double)) ||
    !(x = NAG_ALLOC((n + 1), double))
    )
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Set up solution grid:
   * nag_ode_bvp_ps_lin_cgl_grid (d02ucc).
   * Chebyshev Gauss-Lobatto grid generation.
   */
  nag_ode_bvp_ps_lin_cgl_grid(n, a, b, x, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_ode_bvp_ps_lin_cgl_grid (d02ucc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Set up problem right hand sides for grid. */
  for (i = 0; i < n + 1; i++) f[i] = exact(x[i]);

  scale = 0.5 * (b - a);

  /* Solve on equally spaced grid:
   * nag_ode_bvp_ps_lin_quad_weights (d02uyc).
   * Clenshaw-Curtis quadrature weights for integration using computed
   * Chebyshev coefficients.
   */
  nag_ode_bvp_ps_lin_quad_weights(n, w, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_ode_bvp_ps_lin_quad_weights (d02uyc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Apply the weights, w, to the function values, f, and scale. */
  integ = 0.0;
  for (i = 0; i < n+1; i++) integ = integ + w[i]*f[i];
  integ = scale*integ;

  /* Print function values and weights if required. */
  if (reqwgt) {
    printf("f(x) and integral weights\n\n");
    printf("       x         f(x)        w\n");
    for (i = 0; i < n + 1; i++) {
      printf("%10.4f %10.4f %10.4f\n", x[i], f[i], w[i]);
    }
    printf("\n");
  }

  /* Print approximation to integral. */
  printf("Integral of f(x) from %6.1f to %6.1f = %13.5f\n", a, b, integ);
  if (reqerr) {
    uerr = fabs(integ - 28.0);
    printf("Integral is within a multiple ");
    printf("%8"NAG_IFMT " ", 10 * ((Integer) (uerr/teneps) + 1));
    printf(" of machine precision.\n");
  }
 END:
  if (f) NAG_FREE(f);
  if (w) NAG_FREE(w);
  if (x) NAG_FREE(x);
  return exit_status;
}

static double exact(double x)
{
  return 3.0 * pow(x, 2);
}