/* nag_2d_scat_interpolant(e01sac) Example Program
 *
 * Copyright 1996 Numerical Algorithms Group.
 *
 * Mark 4, 1996.
 * Mark 8 revised, 2004.
 * *********   This example illustrates how to interpolate using Shepherd's method only *************
 */

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

int main(void)
{
  
  Integer exit_status=0, i, j, m, n, nx, ny;
  Integer *triang=0;
  double *grads=0;
  NagError fail;
  Nag_2d_Scat_Method method;
  Nag_Scat_Struct comm;
  double *pf=0, *px=0, *py=0, xhi, xlo, yhi, ylo;
  double x[] = {11.16, 12.85, 19.85, 19.72, 15.91, 0.00, 20.87, 3.45, 14.26, 17.43,
		22.80, 7.58, 25.00, 0.00, 9.66, 5.22, 17.25, 25.00, 12.13, 22.23,
		11.52, 15.20, 7.54, 17.32, 2.14, 0.51, 22.69, 5.47, 21.67, 3.31};
  double y[] = {1.24, 3.06, 10.72, 1.39, 7.74, 20.00, 20.00, 12.78, 17.87, 3.46,
		12.39, 1.98, 11.87, 0.00, 20.00, 14.66, 19.57, 3.87, 10.79, 6.21,
		8.53, 0.00, 10.69, 13.78, 15.03, 8.37, 19.63, 17.13, 14.36, 0.33};
  
  double f[] = {22.15, 22.11, 7.97, 16.83, 15.30, 34.60, 5.74, 41.24, 10.74,
		18.60, 5.47, 29.87, 4.40, 58.20, 4.73, 40.36, 6.43, 8.74,
		13.71, 10.25, 15.74, 21.60, 19.31, 12.11, 53.10, 49.43,
		3.25, 28.63, 5.52, 44.08};
  INIT_FAIL(fail);
  Vprintf("e01sac Example Program Results\n");
  
  /* The number of nodes. */
  m = 30;
      if ( !( triang = NAG_ALLOC(7*m, Integer)) ||
	   !( grads = NAG_ALLOC(2*m, double)) ) 
	{
	  Vprintf("Allocation failure\n");
	  exit_status = -1;
	  goto END;
	}
  
  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.
   */
  
#ifdef E01SAC
  e01sac(method, m, x, y, f, &comm, (Nag_E01_Opt *)0,
	 &fail);
#else
  e01sjc(m, x, y, f, triang, grads, &fail);
#endif
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from e01sac/e01sjc.\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  
  /* Input the domain for evaluating the interpolant. */
  nx = 7;
  xlo = 3.0;
  xhi = 21.0;
  ny = 6;
  ylo = 2.0;
  yhi =  17.0;
  
  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).
   */
  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;
	}
    }
#ifdef E01SAC  
  e01sbc(&comm, nx*ny, px, py, pf, &fail);
#else
  for (i=0; i<nx*ny;i++)
    {
      e01skc(m, x, y, f, triang, grads, px[i], py[i], &pf[i], &fail);
    }
#endif
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from e01sbc/e01skc.\n%s\n", fail.message);
      exit_status = 1;
	  goto END;
    }
  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. */
#ifdef E01SAC  
  e01szc(&comm);
#endif
  if (pf) NAG_FREE(pf);
  if (px) NAG_FREE(px);
  if (py) NAG_FREE(py);
  return exit_status;
}