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

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

int main(int argc, char *argv[])
{
  FILE               *fpin, *fpout;

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

  INIT_FAIL(fail);

  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);
  fprintf(fpout, "nag_2d_scat_interpolant (e01sac) Example Program Results\n");
  /* Skip heading in data file */
  fscanf(fpin, "%*[^\n]");
  /* Input the number of nodes. */
  fscanf(fpin, "%ld", &m);

  if (m >= 3)
    {
      if (!(f = NAG_ALLOC(m, double)) ||
          !(x = NAG_ALLOC(m, double)) ||
          !(y = NAG_ALLOC(m, double)))
        {
          fprintf(fpout, "Allocation failure\n");
          exit_status = -1;
          goto ENDL;
        }
    }
  else
    {
      fprintf(fpout, "Invalid m.\n");
      exit_status = 1;
      return exit_status;
    }
  /* Input the nodes (x,y) and heights, f. */
  for (i = 0; i < m; ++i)
    fscanf(fpin, "%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)
        {
          fprintf(fpout, "\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)
            {
              fprintf(fpout,
                      "Error from nag_2d_scat_interpolant (e01sac).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }
        }
      else if (method == Nag_Shep)
        {
          fprintf(fpout, "\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)
            {
              fprintf(fpout,
                      "Error from nag_2d_scat_interpolant (e01sac).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }

          fprintf(fpout, "   optional.rnw =%8.2f  optional.rnq =%8.2f\n\n",
                  optional.rnw, optional.rnq);
          fprintf(fpout, "   minimum number of data points that lie within "
                  "radius optional.rnq =%3ld\n", optional.minnq);
        }
      /* Input the domain for evaluating the interpolant. */
      fscanf(fpin, "%ld%lf%lf", &nx, &xlo, &xhi);
      fscanf(fpin, "%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)))
            {
              fprintf(fpout, "Allocation failure\n");
              exit_status = -1;
              goto END;
            }
        }
      else
        {
          fprintf(fpout, "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)
            {
              fprintf(fpout, "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);
        }
      fprintf(fpout, "\n          x");
      for (i = 0; i < nx; i++)
        fprintf(fpout, "%8.2f", px[i]);
      fprintf(fpout, "\n     y\n");
      for (i = ny-1; i >= 0; --i)
        {
          fprintf(fpout, "%8.2f   ", py[nx * i]);
          for (j = 0; j < nx; j++)
            fprintf(fpout, "%8.2f", pf[nx * i + j]);
          fprintf(fpout, "\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 (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  if (f) NAG_FREE(f);
  if (x) NAG_FREE(x);
  if (y) NAG_FREE(y);
  return exit_status;
}