/* nag_2d_scat_interpolant(e01sac) Example Program
 *
 * Copyright 1996 Numerical Algorithms Group.
 *
 * Mark 4, 1996.
 */

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

#define MMAX  100
#define NMAX  10
#define NXYM  NMAX*NMAX
#define TDX   MMAX
#define TDQ   MMAX+1

int main(void)
{

  Integer i, isel, j, m, n, nx, ny;
  double  xhi, xlo, yhi, ylo;
  double  f[MMAX], pf[NXYM], px[NXYM], py[NXYM], x[MMAX], y[MMAX];
  Nag_Scat_Struct  comm;
  Nag_2d_Scat_Method method;
  Nag_E01_Opt optional;

  Vprintf("e01sac Example Program Results\n");
  /* Skip heading in data file */
  Vscanf("%*[^\n]");
  /* Input the number of nodes. */
  Vscanf("%ld", &m);

  if (m>=1 && m<=MMAX)
    {
      /* Input the nodes (x,y) and heights, f. */
      for (i=0; i<m; ++i)
        Vscanf("%lf%lf%lf", &x[i], &y[i], &f[i]);

      /* Input the domain for evaluating the interpolant. */
      Vscanf("%ld%lf%lf", &nx, &xlo, &xhi);
      Vscanf("%ld%lf%lf", &ny, &ylo, &yhi);

      for (isel=1; isel<=2; ++isel)
        {
          /* Select the method of interpolation. */
          if (isel==1)
            method = Nag_RC;
          else
            method = Nag_Shep;

          if (method==Nag_RC)
            {
              Vprintf("\nExample 1: Surface interpolation by method "
                      "of Renka and Cline.\n\n");
              /*
               * Generate the triangulation and gradients using the selected
               * method.
               */

              e01sac(method, m, x, y, f, &comm, (Nag_E01_Opt *)0,
                     NAGERR_DEFAULT);
            }
          else if (method==Nag_Shep)
            {
              Vprintf("\n\nExample 2: Surface interpolation by modified "
                      "Shepard's method.\n\n");
              /* Compute the nodal function coefficients. */
              optional.nq = 24;
              optional.nw = 12;
              optional.rnq = -1.0;

              e01sac(method, m, x, y, f, &comm, &optional, NAGERR_DEFAULT);

              Vprintf("   optional.rnw =%8.2f  optional.rnq =%8.2f\n\n",
                      optional.rnw, optional.rnq);
              Vprintf("   minimum number of data points that lie within "
                      "radius optional.rnq =%3ld\n", optional.minnq);
            }
          /*
           * Evaluate the interpolant on a rectangular grid at nx*ny points
           * over the domain (xlo to xhi) x (ylo to yhi).
           */
          n = 0;
          if (nx*ny<=NXYM)
            {
              for (j=0; j<ny; ++j)
                {
                  for (i=0; i<nx; ++i)
                    {
                      px[i+nx*j] = ((double)(nx-i-1) / (nx-1)) * xlo +
                        ((double)(i) / (nx-1)) * xhi;
                      py[i+nx*j] = ((double)(ny-j-1) / (ny-1)) * ylo +
                        ((double)(j) / (ny-1)) * yhi;
                      ++n;
                    }
                }
              if (method == Nag_RC)
                e01sbc(&comm, n, px, py, pf, NAGERR_DEFAULT);
              else if (method == Nag_Shep)
                e01sbc(&comm, n, px, py, pf, NAGERR_DEFAULT);

              Vprintf("\n          x");
              for (i = 0; i < nx; i++)
                Vprintf("%8.2f", px[i]);
              Vprintf("\n     y\n");
              for (i = ny-1; i >= 0; --i)
                {
                  Vprintf("%8.2f   ", py[nx * i]);
                  for (j = 0; j < nx; j++)
                    Vprintf("%8.2f", pf[nx * i + j]);
                  Vprintf("\n");
                }
            }
          /* Free the memory allocated to the pointers in structure comm. */
          e01szc(&comm);
        }
      return EXIT_SUCCESS;
    }
  else
    {
      Vfprintf(stderr, "m is out of range: m = %-3ld\n", m);
      return EXIT_FAILURE;
    }
}