/* nag_2d_triang_interp (e01sjc) Example Program.
 *
 * Copyright 2005 Numerical Algorithms Group.
 *
 * Mark 8, 2004.
 */

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

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

  /* Scalars */
  double   xhi, xlo, yhi, ylo;
  Integer  exit_status, i, j, m, nx, ny;

  /* Arrays */
  double   *f = 0, *grads = 0, *pf = 0, *px = 0, *py = 0, *x = 0, *y = 0;
  Integer  *triang = 0;

  /* Nag Types */
  NagError fail;

  exit_status = 0;
  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_triang_interp (e01sjc) Example Program Results\n");

  /* Skip heading in data file */
  fscanf(fpin, "%*[^\n] ");

  /* Input the number of nodes. */
  fscanf(fpin, "%ld%*[^\n] ", &m);
  if (m >= 3)
    {
      /* Allocate memory */
      if (!(f = NAG_ALLOC(m, double)) ||
          !(grads = NAG_ALLOC(2 * m, double)) ||
          !(x = NAG_ALLOC(m, double)) ||
          !(y = NAG_ALLOC(m, double)) ||
          !(triang = NAG_ALLOC(7*m, Integer)))
        {
          fprintf(fpout, "Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  else
    {
      fprintf(fpout, "Invalid m.\n");
      exit_status = 1;
      goto END;
    }
  /* Input the nodes (X,Y) and heights, F. */
  for (i = 1; i <= m; ++i)
    {
      fscanf(fpin, "%lf%lf%lf%*[^\n] ", &x[i - 1], &y[i - 1], &f[i - 1]);
    }
  /* Generate the triangulation and gradients. */

  /* nag_2d_triang_interp (e01sjc).
   * 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_triang_interp(m, x, y, f, triang, grads, &fail);
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_2d_triang_interp (e01sjc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  /* Evaluate the interpolant on a rectangular grid at NX*NY points */
  /* over the domain (XLO to XHI) x (YLO to YHI). */

  fscanf(fpin, "%ld%lf%lf%*[^\n] ", &nx, &xlo, &xhi);
  fscanf(fpin, "%ld%lf%lf%*[^\n] ", &ny, &ylo, &yhi);

  if (nx > 0 && ny > 0)
    {
      /* Allocate memory */
      if (!(pf = NAG_ALLOC(nx, double)) ||
          !(px = NAG_ALLOC(nx, double)) ||
          !(py = NAG_ALLOC(ny, double)))
        {
          fprintf(fpout, "Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  else
    {
      fprintf(fpout, "Invalid nx or ny.\n");
      exit_status = 1;
      goto END;
    }

  for (i = 1; i <= nx; ++i)
    {
      px[i - 1] = (double)(nx - i) / (nx - 1) * xlo +
                  (double)(i - 1) / (nx - 1) * xhi;
    }
  for (i = 1; i <= ny; ++i)
    {
      py[i - 1] = (double)(ny - i) / (ny - 1) * ylo +
                  (double)(i - 1) / (ny - 1) * yhi;
    }
  fprintf(fpout, "\n");
  fprintf(fpout, "%s", "          X");
  for (i = 1; i <= nx; ++i)
    {
      fprintf(fpout, "%8.2f", px[i - 1]);
    }
  fprintf(fpout, "\n");
  fprintf(fpout, "%s", "     Y");
  fprintf(fpout, "\n");
  for (i = ny; i >= 1; --i)
    {
      for (j = 1; j <= nx; ++j)
        {
          /* nag_2d_triang_eval (e01skc).
           * A function to evaluate, at a set of points, the
           * two-dimensional interpolant function generated by
           * nag_2d_scat_interpolant (e01sac)
           */
          nag_2d_triang_eval(m, x, y, f, triang, grads, px[j - 1],
                             py[i - 1], &pf[j - 1], &fail);
          if (fail.code != NE_NOERROR)
            {
              fprintf(fpout, "Error from nag_2d_triang_eval (e01skc).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }
        }
      fprintf(fpout, "%8.2f", py[i - 1]);
      fprintf(fpout, "%3s", "");
      for (j = 1; j <= nx; ++j)
        {
          fprintf(fpout, "%8.2f", pf[j - 1]);
        }
      fprintf(fpout, "\n");
    }
 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  if (f) NAG_FREE(f);
  if (grads) NAG_FREE(grads);
  if (pf) NAG_FREE(pf);
  if (px) NAG_FREE(px);
  if (py) NAG_FREE(py);
  if (x) NAG_FREE(x);
  if (y) NAG_FREE(y);
  if (triang) NAG_FREE(triang);

  return exit_status;
}