/* nag_2d_triang_interp (e01sjc) Example Program.

 *

 * Copyright 2005 Numerical Algorithms Group.

 *

 * Mark 8, 2004.

 */



#include <stdio.h>

#include <nag.h>

#include <nag_stdlib.h>

#include <nage01.h>



int main(void)

{



  /* 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);

  Vprintf( "nag_2d_triang_interp (e01sjc) Example Program Results\n");



  /* Skip heading in data file */

  Vscanf("%*[^\n] ");



  /* Input the number of nodes. */

  Vscanf("%ld%*[^\n] ",   &m);

  if (m >= 3)

    {

      /* Allocate memory */

      if ( !(f = NAG_ALLOC(m, double)) ||

           !(grads = NAG_ALLOC(2 * m, double)) ||

           !(pf = NAG_ALLOC(25, double)) ||

           !(px = NAG_ALLOC(25, double)) ||

           !(py = NAG_ALLOC(25, double)) ||

           !(x = NAG_ALLOC(m, double)) ||

           !(y = NAG_ALLOC(m, double)) ||

           !(triang = NAG_ALLOC(7*m, Integer)) )

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }



    }

  else

    {

      Vprintf("Invalid m.\n");

      exit_status = 1;

      return exit_status;

    }

  /* Input the nodes (X,Y) and heights, F. */

  for (i = 1; i <= m; ++i)

    {

      Vscanf("%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)

    {

      Vprintf("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). */



  Vscanf("%ld%lf%lf%*[^\n] ",   &nx,  &xlo,  &xhi);

  Vscanf("%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))  )

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }

    }

  else

    {

      Vprintf("Invalid nx or ny.\n");

      exit_status = 1;

      return exit_status;

    }



  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;

    }

  Vprintf("\n");

  Vprintf("%s", "          X");

  for (i = 1; i <= nx; ++i)

    {

      Vprintf("%8.2f",   px[i - 1]);

    }

  Vprintf("\n");

  Vprintf("%s", "     Y");

  Vprintf("\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)

            {

              Vprintf("Error from nag_2d_triang_eval (e01skc).\n%s\n",

                      fail.message);

              exit_status = 1;

              goto END;

            }

        }

      Vprintf("%8.2f", py[i - 1]);

      Vprintf("%3s", "");

      for (j = 1; j <= nx; ++j)

        {

          Vprintf("%8.2f",   pf[j - 1]);

        }

      Vprintf("\n");

    }

 END:

  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;

}