/* nag_inteq_fredholm2_smooth (d05abc) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.5, 2022.
 */
#include <math.h>
#include <nag.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 */
  static double ruser[2] = {-1.0, -1.0};
  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");

  /* For communication with user-supplied functions: */
  comm.user = ruser;

  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_math_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 = %" NAG_IFMT "\n\n", n);
    printf("Solution on first %2" NAG_IFMT " 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_math_pi * (double)(i) / (double)(2 * n - 1));
      printf("%3" NAG_IFMT "%15.5f%15.5f%15.5e\n", i, y, f[i], c[i]);
    }
    printf("\n");

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

    if (fail.code != NE_NOERROR) {
      printf("Error from nag_sum_chebyshev (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");

    NAG_FREE(c);
    NAG_FREE(f);
  }

END:

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

  return exit_status;
}

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

  if (comm->user[0] == -1.0) {
    printf("(User-supplied callback k, first invocation.)\n");
    comm->user[0] = 0.0;
  }
  return alpha / (pow(alpha, 2) + pow(x - s, 2));
}

static double NAG_CALL g(double x, Nag_Comm *comm) {
  if (comm->user[1] == -1.0) {
    printf("(User-supplied callback g, first invocation.)\n");
    comm->user[1] = 0.0;
  }
  return 1.0;
}