/* nag_2d_spline_interpolant (e01dac) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.2, 2017.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nage01.h>
#include <nage02.h>

#define F(I, J)  f[my*(I)+(J)]
#define FG(I, J) fg[npy*(I)+(J)]
#define C(I, J)  spline.c[my*(I)+(J)]
int main(void)
{
  Integer exit_status = 0, i, j, mx, my, npx, npy;
  NagError fail;
  Nag_2dSpline spline;
  double *f = 0, *fg = 0, step, *tx = 0, *ty = 0, *x = 0, xhi, xlo;
  double *y = 0, yhi, ylo;

  INIT_FAIL(fail);

  /* Initialize spline */
  spline.lamda = 0;
  spline.mu = 0;
  spline.c = 0;

  printf("nag_2d_spline_interpolant (e01dac) Example Program Results\n");
  scanf("%*[^\n]"); /* Skip heading in data file */
  /* Read the number of x points, mx, and the values of the
   * x co-ordinates.
   */
  scanf("%" NAG_IFMT "%" NAG_IFMT "", &mx, &my);
  if (mx >= 4 && my >= 4) {
    if (!(f = NAG_ALLOC(mx * my, double)) ||
        !(x = NAG_ALLOC(mx, double)) || !(y = NAG_ALLOC(my, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  }
  else {
    printf("Invalid mx or my.\n");
    exit_status = 1;
    return exit_status;
  }
  for (i = 0; i < mx; i++)
    scanf("%lf", &x[i]);
  /* Read the number of y points, my, and the values of the
   * y co-ordinates.
   */
  for (i = 0; i < my; i++)
    scanf("%lf", &y[i]);
  /* Read the function values at the grid points. */
  for (j = 0; j < my; j++)
    for (i = 0; i < mx; i++)
      scanf("%lf", &F(i, j));
  /* Generate the (x,y,f) interpolating bicubic B-spline. */
  /* nag_2d_spline_interpolant (e01dac).
   * Interpolating function, bicubic spline interpolant, two
   * variables
   */
  nag_2d_spline_interpolant(mx, my, x, y, f, &spline, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_2d_spline_interpolant (e01dac).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Print the knot sets, lamda and mu. */
  printf("Distinct knots in x direction located  at\n");
  for (j = 3; j < spline.nx - 3; j++)
    printf("%12.4f%s", spline.lamda[j],
           ((j - 3) % 5 == 4 || j == spline.nx - 4) ? "\n" : " ");
  printf("\nDistinct knots in y direction located  at\n");
  for (j = 3; j < spline.ny - 3; j++)
    printf("%12.4f%s", spline.mu[j],
           ((j - 3) % 5 == 4 || j == spline.ny - 4) ? "\n" : " ");
  /* Print the spline coefficients. */
  printf("\nThe B-Spline coefficients:\n");
  for (i = 0; i < mx; i++) {
    for (j = 0; j < my; j++)
      printf("%9.4f", C(i, j));
    printf("\n");
  }

  /* Evaluate the spline on a regular rectangular grid at npx*npy
   * points over the domain (xlo to xhi) x (ylo to yhi).
   */
  scanf("%" NAG_IFMT "%lf%lf", &npx, &xlo, &xhi);
  scanf("%" NAG_IFMT "%lf%lf", &npy, &ylo, &yhi);
  if (npx >= 1 && npy >= 1) {
    if (!(fg = NAG_ALLOC(npx * npy, double)) ||
        !(tx = NAG_ALLOC(npx, double)) || !(ty = NAG_ALLOC(npy, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  }
  else {
    printf("Invalid npx or npy.\n");
    exit_status = 1;
    return exit_status;
  }
  step = (xhi - xlo) / (double) (npx - 1);
  printf("\nSpline evaluated on a regular mesh "
         "     (x across, y down): \n     ");
  /* Generate nx equispaced x co-ordinates. */
  for (i = 0; i < npx; i++) {
    tx[i] = MIN(xlo + i * step, xhi);
    printf("   %5.2f ", tx[i]);
  }
  step = (yhi - ylo) / (npy - 1);
  for (i = 0; i < npy; i++)
    ty[i] = MIN(ylo + i * step, yhi);

  /* Evaluate the spline. */
  /* nag_2d_spline_eval_rect (e02dfc).
   * Evaluation of bicubic spline, at a mesh of points
   */
  nag_2d_spline_eval_rect(npx, npy, tx, ty, fg, &spline, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_2d_spline_eval_rect (e02dfc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Print the results. */
  printf("\n");
  for (j = 0; j < npy; j++) {
    printf("%5.2f", ty[j]);
    for (i = 0; i < npx; i++)
      printf("%8.3f ", FG(i, j));
    printf("\n");
  }
  /* Free memory allocated by nag_2d_spline_interpolant (e01dac) */
END:
  NAG_FREE(spline.lamda);
  NAG_FREE(spline.mu);
  NAG_FREE(spline.c);
  NAG_FREE(f);
  NAG_FREE(x);
  NAG_FREE(y);
  NAG_FREE(fg);
  NAG_FREE(tx);
  NAG_FREE(ty);
  return exit_status;
}