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

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

int main(void)
{

  Integer exit_status=0, i, isel, j, m, n, nx, ny;
  NagError fail;
  Nag_2d_Scat_Method method;
  Nag_E01_Opt optional;
  Nag_Scat_Struct comm;
  double *f=0, *pf=0, *px=0, *py=0, *x=0, xhi, xlo, *y=0, yhi, ylo;

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

  if (m>=3)
    {
      if ( !( f = NAG_ALLOC(m, double)) ||
           !( x = NAG_ALLOC(m, double)) ||
           !( y = NAG_ALLOC(m, double)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto ENDL;
        }
    }
  else
    {
      Vprintf("Invalid m.\n");
      exit_status = 1;
      return exit_status;
    }
  /* Input the nodes (x,y) and heights, f. */
  for (i=0; i<m; ++i)
    Vscanf("%lf%lf%lf", &x[i], &y[i], &f[i]);

  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.
           */

          /* nag_2d_scat_interpolant (e01sac).
           * 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_2d_scat_interpolant(method, m, x, y, f, &comm, (Nag_E01_Opt *)0,
                                  &fail);
          if (fail.code != NE_NOERROR)
            {
              Vprintf("Error from nag_2d_scat_interpolant (e01sac).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }
        }
      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;

          /* nag_2d_scat_interpolant (e01sac), see above. */
          nag_2d_scat_interpolant(method, m, x, y, f, &comm, &optional, &fail);
          if (fail.code != NE_NOERROR)
            {
              Vprintf("Error from nag_2d_scat_interpolant (e01sac).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }

          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);
        }
      /* Input the domain for evaluating the interpolant. */
      Vscanf("%ld%lf%lf", &nx, &xlo, &xhi);
      Vscanf("%ld%lf%lf", &ny, &ylo, &yhi);

      if (nx>=1 && ny>=1)
        {
          if ( !( pf = NAG_ALLOC(nx*ny, double)) ||
               !( px = NAG_ALLOC(nx*ny, double)) ||
               !( py = NAG_ALLOC(nx*ny, double)) )
            {
              Vprintf("Allocation failure\n");
              exit_status = -1;
              goto END;
            }
        }
      else
        {
          Vprintf("Invalid nx or ny.\n");
          exit_status = 1;
          return exit_status;
        }
      /*
       * Evaluate the interpolant on a rectangular grid at nx*ny points
       * over the domain (xlo to xhi) x (ylo to yhi).
       */
      n = 0;
      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)
        {
          /* nag_2d_scat_eval (e01sbc).
           * A function to evaluate, at a set of points, the
           * two-dimensional interpolant function generated by
           * nag_2d_scat_interpolant (e01sac)
           */
          nag_2d_scat_eval(&comm, n, px, py, pf, &fail);
          if (fail.code != NE_NOERROR)
            {
              Vprintf("Error from nag_2d_scat_eval (e01sbc).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }
        }
      else if (method == Nag_Shep)
        {
          /* nag_2d_scat_eval (e01sbc), see above. */
          nag_2d_scat_eval(&comm, n, px, py, pf, &fail);
        }
      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");
        }
    END:
      /* Free the memory allocated to the pointers in structure comm. */
      /* nag_2d_scat_free (e01szc).
       * Freeing function for use with nag_2d_scat_eval (e01sbc)
       */
      nag_2d_scat_free(&comm);
      if (pf) NAG_FREE(pf);
      if (px) NAG_FREE(px);
      if (py) NAG_FREE(py);
    }
 ENDL:
  if (f) NAG_FREE(f);
  if (x) NAG_FREE(x);
  if (y) NAG_FREE(y);
  return exit_status;
}