/* nag_inteq_fredholm2_smooth (d05abc) Example Program.
 *
 * Copyright 2011, Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagc06.h>
#include <nagd05.h>
#include <nagx01.h>

#ifdef __cplusplus
extern "C" {
#endif
  static double NAG_CALL k(double x, double s, Nag_Comm *comm);
  static double NAG_CALL g(double x, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{
  /* Scalars */
  double      a = -1.0, b = 1.0;
  double      lambda, x0;
  Integer     exit_status = 0;
  Integer     i, lx, n;
  Nag_Boolean ev = Nag_TRUE, odorev = Nag_TRUE;
  /* Arrays */
  double      *c = 0, *chebr = 0, *f = 0, *x = 0;
  /* NAG types */
  Nag_Comm comm;
  NagError fail;
  Nag_Series s = Nag_SeriesEven;

  INIT_FAIL(fail);

  printf("nag_inteq_fredholm2_smooth (d05abc) Example Program Results\n");

  x0 = 0.5 * (a + b);

  /* Set up uniform grid to evaluate Chebyshev polynomials. */
  lx = (Integer) (4.000001 * (b - x0)) + 1;

  if (
      !(x = NAG_ALLOC(lx, double)) ||
      !(chebr = NAG_ALLOC(lx, double))
      )
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  x[0] = x0;
  for (i = 1; i < lx; i++)
    x[i] = x[i - 1] + 0.25;

  printf("\nSolution is even\n");

  lambda = -1.0 / nag_pi;

  for (n = 5; n <= 10; n += 5)
    {
      if (
          !(f = NAG_ALLOC(n, double)) ||
          !(c = NAG_ALLOC(n, double))
          )
        {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }

      /*
        nag_inteq_fredholm2_smooth (d05abc).
        Linear non-singular Fredholm integral equation, second kind,
        smooth kernel.
      */
      nag_inteq_fredholm2_smooth(lambda, a, b, n, k, g, odorev, ev, f, c, &comm,
                                 &fail);

      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_inteq_fredholm2_smooth (d05abc).\n%s\n",
                 fail.message);
          exit_status = 1;
          goto END;
        }

      printf("\nResults for n = %ld\n\n", n);
      printf("Solution on first %2ld Chebyshev points and Chebyshev"
	     " coefficients\n",n);
      printf("%3s%12s%18s%12s\n","i","x","f[i]","c[i]");
      for (i = 0; i < n; i++)
        {
          double y = cos(nag_pi*(double) (i)/(double)(2*n-1));
          printf("%3ld%15.5f%15.5f%15.5e\n", i, y, f[i], c[i]);
        }
      printf("\n");

      /*
        Evaluate and print solution on uniform grid.
        nag_sum_cheby_series (c06dcc).
        Sum of a Chebyshev series at a set of points.
      */
      nag_sum_cheby_series(x, lx, a, b, c, n, s, chebr, &fail);

      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_sum_cheby_series (c06dcc).\n%s\n",
                 fail.message);
          exit_status = 1;
          goto END;
        }

      printf("Solution on evenly spaced grid\n");
      printf("\n     x           f(x)\n");
      for (i = 0; i < lx; i++)
        printf("%8.4f%15.5f\n", x[i], chebr[i]);
      printf("\n");

      if (c) NAG_FREE(c);
      if (f) NAG_FREE(f);
    }

 END:

  if (c) NAG_FREE(c);
  if (f) NAG_FREE(f);
  if (chebr) NAG_FREE(chebr);
  if (x) NAG_FREE(x);

  return exit_status;
}

static double NAG_CALL k(double x, double s, Nag_Comm *comm)
{
  /* Scalars */
  double alpha = 1.0;

  return alpha / (pow(alpha, 2) + pow(x - s, 2));
}
static double NAG_CALL g(double x, Nag_Comm *comm)
{
  return 1.0;
}