/* nag_interp_dim2_triang_interp (e01sjc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.0, 2023.
 */

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

int main(void) {

  /* Scalars */
  double xhi, xlo, yhi, ylo;
  Integer exit_status, i, j, m, nx, ny;

  /* Arrays */
  double *f = 0, *grads = 0, *pf = 0, *px = 0, *py = 0, *x = 0, *y = 0;
  Integer *triang = 0;

  /* Nag Types */
  NagError fail;

  exit_status = 0;
  INIT_FAIL(fail);

  printf("nag_interp_dim2_triang_interp (e01sjc) Example Program Results\n");

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

  /* Input the number of nodes. */
  scanf("%" NAG_IFMT "%*[^\n] ", &m);
  if (m >= 3) {
    /* Allocate memory */
    if (!(f = NAG_ALLOC(m, double)) || !(grads = NAG_ALLOC(2 * m, double)) ||
        !(x = NAG_ALLOC(m, double)) || !(y = NAG_ALLOC(m, double)) ||
        !(triang = NAG_ALLOC(7 * m, Integer))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  } else {
    printf("Invalid m.\n");
    exit_status = 1;
    goto END;
  }
  /* Input the nodes (X,Y) and heights, F. */
  for (i = 1; i <= m; ++i) {
    scanf("%lf%lf%lf%*[^\n] ", &x[i - 1], &y[i - 1], &f[i - 1]);
  }
  /* Generate the triangulation and gradients. */

  /* nag_interp_dim2_triang_interp (e01sjc).
   * A function to generate a two-dimensional surface
   * interpolating a set of data points, using either the
   * method of Renka and Cline or the modified Shepard's
   * method
   */
  nag_interp_dim2_triang_interp(m, x, y, f, triang, grads, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_interp_dim2_triang_interp (e01sjc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Evaluate the interpolant on a rectangular grid at NX*NY points */
  /* over the domain (XLO to XHI) x (YLO to YHI). */
  scanf("%" NAG_IFMT "%lf%lf%*[^\n] ", &nx, &xlo, &xhi);
  scanf("%" NAG_IFMT "%lf%lf%*[^\n] ", &ny, &ylo, &yhi);

  if (nx > 0 && ny > 0) {
    /* Allocate memory */
    if (!(pf = NAG_ALLOC(nx, double)) || !(px = NAG_ALLOC(nx, double)) ||
        !(py = NAG_ALLOC(ny, double))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  } else {
    printf("Invalid nx or ny.\n");
    exit_status = 1;
    goto END;
  }

  for (i = 1; i <= nx; ++i) {
    px[i - 1] =
        (double)(nx - i) / (nx - 1) * xlo + (double)(i - 1) / (nx - 1) * xhi;
  }
  for (i = 1; i <= ny; ++i) {
    py[i - 1] =
        (double)(ny - i) / (ny - 1) * ylo + (double)(i - 1) / (ny - 1) * yhi;
  }
  printf("\n");
  printf("%s", "          X");
  for (i = 1; i <= nx; ++i) {
    printf("%8.2f", px[i - 1]);
  }
  printf("\n");
  printf("%s", "     Y");
  printf("\n");
  for (i = ny; i >= 1; --i) {
    for (j = 1; j <= nx; ++j) {
      /* nag_interp_dim2_triang_eval (e01skc).
       * A function to evaluate, at a set of points, the
       * two-dimensional interpolant function generated by
       * nag_interp_dim2_triang_interp (e01sjc).
       */
      nag_interp_dim2_triang_eval(m, x, y, f, triang, grads, px[j - 1],
                                  py[i - 1], &pf[j - 1], &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_interp_dim2_triang_eval (e01skc).\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }
    }
    printf("%8.2f", py[i - 1]);
    printf("%3s", "");
    for (j = 1; j <= nx; ++j) {
      printf("%8.2f", pf[j - 1]);
    }
    printf("\n");
  }
END:
  NAG_FREE(f);
  NAG_FREE(grads);
  NAG_FREE(pf);
  NAG_FREE(px);
  NAG_FREE(py);
  NAG_FREE(x);
  NAG_FREE(y);
  NAG_FREE(triang);

  return exit_status;
}