/* 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 =%3"NAG_IFMT"\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;

}