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

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

#define MXMAX 20
#define MYMAX 20
#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 i, j, mx, my, npx, npy;
  double  f[MXMAX*MYMAX], x[MXMAX], y[MYMAX];
  double fg[MXMAX*MYMAX], tx[MXMAX], ty[MYMAX];
  double  xhi, yhi, xlo, ylo, step;
  Nag_2dSpline spline;

  Vprintf("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>MXMAX ||  my>MYMAX)
    {
      Vfprintf(stderr, "mx or my is out of range: mx = %5ld\n,"
               "my = %5ld\n",mx,my);
      return EXIT_FAILURE;
    }
  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. */
  e01dac(mx, my, x, y, f, &spline, NAGERR_DEFAULT);

  /* 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<=MXMAX && npy<=MYMAX)
    {
      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. */
      e02dfc(npx, npy, tx, ty, fg, &spline, NAGERR_DEFAULT);

      /* 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 e01dac */
      NAG_FREE(spline.lamda);
      NAG_FREE(spline.mu);
      NAG_FREE(spline.c);
      return EXIT_SUCCESS;
    }
  else
    {
      Vfprintf(stderr, "npx or npy is out of range: npx = %5ld, npy = %5ld\n",
               npx,npy);
      /* Free memory allocated by e01dac */
      NAG_FREE(spline.lamda);
      NAG_FREE(spline.mu);
      NAG_FREE(spline.c);
      return EXIT_FAILURE;
    }
}