```/* 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;
}
```