/* nag_ode_bvp_ps_lin_cheb_eval (d02uzc) Example Program.
 *
 * Copyright 2011, Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */

#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd02.h>
#include <nagx01.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, k, m, n;
  double       a = -0.24 * nag_pi, b = 0.5 * nag_pi;
  double       deven, dmap, fseries, t, uerr, xeven, xmap;
  double       teneps = 10.0 * nag_machine_precision;
  /*  Arrays */
  double      *c = 0, *f = 0, *x = 0;
  /* NAG types */
  Nag_Boolean reqerr = Nag_FALSE;
  NagError    fail;

  INIT_FAIL(fail);

  printf("nag_ode_bvp_ps_lin_cheb_eval (d02uzc) Example Program Results \n\n");

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

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

  /* Set up Chebyshev 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;
  }

  /* Evaluate function on grid and get interpolating Chebyshev coefficients. */
  for (i = 0; i < n + 1; i++) f[i] = exact(x[i]);

  /* nag_ode_bvp_ps_lin_coeffs (d02uac).
   * Coefficients of Chebyshev interpolating polynomial
   * from function values on Chebyshev grid.
   */
  nag_ode_bvp_ps_lin_coeffs(n, f, c, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_ode_bvp_ps_lin_coeffs (d02uac).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Evaluate Chebyshev series manually by evaluating each Chebyshev
   * polynomial in turn at new equispaced (m+1) grid points.
   * Chebyshev series on [-1,1] map of [a,b].
  */
  xmap = -1.0;
  dmap = 2.0/(double) (m - 1);
  xeven = a;
  deven = (b - a)/(double) (m - 1);
  printf("    x_even      x_map      Sum\n");
  uerr = 0.0;
  for (i = 0; i < m; i++) {
    fseries = 0.0;
    for (k = 0; k < n + 1; k++) {
      /* nag_ode_bvp_ps_lin_cheb_eval (d02uzc).
       * Chebyshev polynomial evaluation, T_k(x).
       */
      nag_ode_bvp_ps_lin_cheb_eval(k, xmap, &t, &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_ode_bvp_ps_lin_cheb_eval (d02uzc).\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }

      fseries = fseries + c[k] * t;
    }
    uerr = MAX(uerr, fabs(fseries - exact(xeven)));
    printf("%10.4f %10.4f %10.4f \n", xeven, xmap, fseries);
    xmap = MIN(1.0, xmap + dmap);
    xeven = xeven + deven;
  }

  if (reqerr) {
    printf("\nError in coefficient sum is < ");
    printf("%8"NAG_IFMT " ", 10 * ((Integer) (uerr/teneps) + 1));
    printf(" * machine precision.\n");
  }
 END:
  if (c) NAG_FREE(c);
  if (f) NAG_FREE(f);
  if (x) NAG_FREE(x);
  return exit_status;
}

static double exact(double x)
{
  return x + exp(-x);
}