/* nag_interp_dim4_scat_shep_eval (e01tlc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.0, 2023.
 */

#include <math.h>
#include <nag.h>
#include <stdio.h>

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

#define X(I, J) x[I * 4 + J]
#define XE(I, J) xe[I * 4 + J]

int main(void) {
  /* Scalars */
  Integer exit_status, i, m, n, nq, nw, liq, lrq, lstate, subid;
  Integer lseed = 1;
  double fun;
  Nag_BaseRNG genid;
  NagError fail;
  /* Arrays */
  double *f = 0, *q = 0, *qx = 0, *rq = 0, *xe = 0, *x = 0;
  Integer *iq = 0, *state = 0;
  Integer seed[1], seed2[1];

  exit_status = 0;

  INIT_FAIL(fail);

  printf("nag_interp_dim4_scat_shep_eval (e01tlc) Example Program Results\n");

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

  /* Input the seeds. */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &seed[0], &seed2[0]);

  /* Choose the base generator */
  genid = Nag_Basic;
  subid = 0;

  /* Get the length of the state array */
  lstate = -1;
  nag_rand_init_repeat(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_init_repeat (g05kfc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Input the number of nodes. */
  scanf("%" NAG_IFMT "%*[^\n] ", &m);

  /* Allocate memory */
  lrq = 21 * m + 11;
  liq = 2 * m + 1;
  if (!(f = NAG_ALLOC(m, double)) || !(x = NAG_ALLOC(m * 4, double)) ||
      !(rq = NAG_ALLOC(lrq, double)) || !(iq = NAG_ALLOC(liq, Integer)) ||
      !(state = NAG_ALLOC(lstate, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Initialize the generator to a repeatable sequence */
  nag_rand_init_repeat(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_init_repeat (g05kfc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Generate the data points X */
  nag_rand_dist_uniform01(m * 4, state, x, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_dist_uniform01 (g05sac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Evaluate F */
  for (i = 0; i < m; ++i) {
    f[i] = funct(&X(i, 0));
  }

  /* Generate the interpolant. */
  nq = 0;
  nw = 0;

  /* nag_interp_dim4_scat_shep (e01tkc).
   * Interpolating functions, modified Shepard's method, four
   * variables
   */
  nag_interp_dim4_scat_shep(m, x, f, nw, nq, iq, rq, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_interp_dim4_scat_shep (e01tkc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Input the number of evaluation points. */
  scanf("%" NAG_IFMT "%*[^\n] ", &n);

  /* Allocate memory for nag_interp_dim4_scat_shep_eval (e01tlc) */
  if (!(q = NAG_ALLOC(n, double)) || !(qx = NAG_ALLOC(n * 4, double)) ||
      !(xe = NAG_ALLOC(n * 4, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Generate repeatable evaluation points. */
  nag_rand_init_repeat(genid, subid, seed2, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_init_repeat (g05kfc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  nag_rand_dist_uniform01(n * 4, state, xe, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_dist_uniform01 (g05sac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* nag_interp_dim4_scat_shep_eval (e01tlc).
   * Interpolated values, evaluate interpolant and first derivatives
   * computed by nag_interp_dim4_scat_shep (e01tkc).
   */
  fail.print = Nag_TRUE;
  nag_interp_dim4_scat_shep_eval(m, x, f, iq, rq, n, xe, q, qx, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_interp_dim4_scat_shep_eval (e01tlc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\n     i     f(x)      q(x)   |f(x)-q(x)|\n");
  for (i = 0; i < n; ++i) {
    fun = funct(&XE(i, 0));
    printf("%6" NAG_IFMT "%10.4f%10.4f%10.4f\n", i, fun, q[i],
           fabs(fun - q[i]));
  }

END:
  NAG_FREE(f);
  NAG_FREE(q);
  NAG_FREE(qx);
  NAG_FREE(rq);
  NAG_FREE(xe);
  NAG_FREE(x);
  NAG_FREE(iq);
  NAG_FREE(state);

  return exit_status;
}

static double NAG_CALL funct(double x[]) {
  /* Scalars */
  double ret_val;

  ret_val = ((1.25 + cos(5.4 * x[3])) * cos(6.0 * x[0]) * cos(6.0 * x[1])) /
            (6.0 * (1.0 + pow((3.0 * x[2] - 1.0), 2.0)));
  return ret_val;
}