/* 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 liq, lrq, nq, nw;
  NagError fail;
  Nag_2d_Scat_Method method;
  Nag_E01_Opt optional;
  Nag_Scat_Struct comm;
  double *pf=0, *px=0, *py=0, *rq=0, *qx=0, *qy=0, xhi, xlo, yhi, ylo;
  Integer *iq=0;
  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;
  n = 5;
      liq = 2*m+1;
      lrq = 6*m+5;
      /* Allocate memory */
      if ( !(rq = NAG_ALLOC(lrq, double)) ||
	   !(iq = NAG_ALLOC(liq, Integer))
	   )
	{
	  Vprintf("Allocation failure\n");
	  exit_status = -1;
	  goto END;
	}

#ifdef E01SAC
        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;
          
          e01sac(method, m, x, y, f, &comm, &optional, &fail);
          
          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);
#else
  nq = 0;
  nw = 0;
  e01sgc(m, x, y, f, nw, nq, iq, rq, &fail);
#endif
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from e01sac/e01sgc.\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)) ||
	       !(qx = NAG_ALLOC(nx*ny, double)) ||
	       !(qy = 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;
            }
        }
#ifdef E01SAC
      e01sbc(&comm, n, px, py, pf, &fail);
#else
      e01shc(m, x, y, f, iq, rq, n, px, py, pf, qx, qy, &fail);
#endif
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from e01sbc/e01shc.\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 (rq) NAG_FREE(rq);
      if (iq) NAG_FREE(rq);
      if (pf) NAG_FREE(pf);
      if (px) NAG_FREE(px);
      if (py) NAG_FREE(py);
      if (qx) NAG_FREE(qx);
      if (qy) NAG_FREE(qy);
  return exit_status;
}