/* nag_2d_spline_interpolant (e01dac) Example Program.
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 *
 * Mark 6 revised, 2000.
 * Mark 8 revised, 2004.
 */

#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, *y=0, yhi, ylo;

  INIT_FAIL(fail);

  Vprintf("nag_2d_spline_interpolant (e01dac) Example Program Results\n");
  Vscanf("%*[^\n]"); /* Skip heading in data file */
  /* Read the number of x points, mx, and the values of the
   * x co-ordinates.
   */
  Vscanf("%ld%ld",&mx, &my);
  if (mx>=4 && my>= 4)
    {
      if ( !( f = NAG_ALLOC(mx*my, double)) ||
           !( x = NAG_ALLOC(mx, double)) ||
           !( y = NAG_ALLOC(my, double)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  else
    {
      Vprintf("Invalid mx or my.\n");
      exit_status = 1;
      return exit_status;
    }
  for (i=0; i<mx;  i++)
    Vscanf("%lf",&x[i]);
  /* Read the number of y points, my, and the values of the
   * y co-ordinates.
   */
  for (i=0; i<my; i++)
    Vscanf("%lf",&y[i]);
  /* Read the function values at the grid points. */
  for (j=0; j<my; j++)
    for (i=0; i<mx; i++)
      Vscanf("%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)
    {
      Vprintf("Error from nag_2d_spline_interpolant (e01dac).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

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

  /* Evaluate the spline on a regular rectangular grid at npx*npy
   * points over the domain (xlo to xhi) x (ylo to yhi).
   */
  Vscanf("%ld%lf%lf",&npx,&xlo,&xhi);
  Vscanf("%ld%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)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  else
    {
      Vprintf("Invalid npx or npy.\n");
      exit_status = 1;
      return exit_status;
    }
  step = (xhi-xlo)/(double)(npx-1);
  Vprintf("\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);
      Vprintf("   %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)
    {
      Vprintf("Error from nag_2d_spline_eval_rect (e02dfc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  /* Print the results. */
  Vprintf("\n");
  for (j=0; j<npy; j++)
    {
      Vprintf("%5.2f",ty[j]);
      for (i=0; i<npx; i++)
        Vprintf("%8.3f ",FG(i,j));
      Vprintf("\n");

    }
  /* Free memory allocated by nag_2d_spline_interpolant (e01dac) */
 END:
  NAG_FREE(spline.lamda);
  NAG_FREE(spline.mu);
  NAG_FREE(spline.c);
  if (f) NAG_FREE(f);
  if (x) NAG_FREE(x);
  if (y) NAG_FREE(y);
  if (fg) NAG_FREE(fg);
  if (tx) NAG_FREE(tx);
  if (ty) NAG_FREE(ty);
  return exit_status;

}