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