```/* nag_2d_triang_interp (e01sjc) Example Program.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.1, 2017.
*/

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

int main(void)
{

/* Scalars */
double xhi, xlo, yhi, ylo;
Integer exit_status, i, j, m, nx, ny;

/* Arrays */
double *f = 0, *grads = 0, *pf = 0, *px = 0, *py = 0, *x = 0, *y = 0;
Integer *triang = 0;

/* Nag Types */
NagError fail;

exit_status = 0;
INIT_FAIL(fail);

printf("nag_2d_triang_interp (e01sjc) Example Program Results\n");

/* Skip heading in data file */
scanf("%*[^\n] ");

/* Input the number of nodes. */
scanf("%" NAG_IFMT "%*[^\n] ", &m);
if (m >= 3) {
/* Allocate memory */
if (!(f = NAG_ALLOC(m, double)) ||
!(grads = NAG_ALLOC(2 * m, double)) ||
!(x = NAG_ALLOC(m, double)) ||
!(y = NAG_ALLOC(m, double)) || !(triang = NAG_ALLOC(7 * m, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
else {
printf("Invalid m.\n");
exit_status = 1;
goto END;
}
/* Input the nodes (X,Y) and heights, F. */
for (i = 1; i <= m; ++i) {
scanf("%lf%lf%lf%*[^\n] ", &x[i - 1], &y[i - 1], &f[i - 1]);
}
/* Generate the triangulation and gradients. */

/* nag_2d_triang_interp (e01sjc).
* A function to generate a two-dimensional surface
* interpolating a set of data points, using either the
* method of Renka and Cline or the modified Shepard's
* method
*/
nag_2d_triang_interp(m, x, y, f, triang, grads, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_2d_triang_interp (e01sjc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}

/* Evaluate the interpolant on a rectangular grid at NX*NY points */
/* over the domain (XLO to XHI) x (YLO to YHI). */
scanf("%" NAG_IFMT "%lf%lf%*[^\n] ", &nx, &xlo, &xhi);
scanf("%" NAG_IFMT "%lf%lf%*[^\n] ", &ny, &ylo, &yhi);

if (nx > 0 && ny > 0) {
/* Allocate memory */
if (!(pf = NAG_ALLOC(nx, double)) ||
!(px = NAG_ALLOC(nx, double)) || !(py = NAG_ALLOC(ny, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
else {
printf("Invalid nx or ny.\n");
exit_status = 1;
goto END;
}

for (i = 1; i <= nx; ++i) {
px[i - 1] = (double) (nx - i) / (nx - 1) * xlo +
(double) (i - 1) / (nx - 1) * xhi;
}
for (i = 1; i <= ny; ++i) {
py[i - 1] = (double) (ny - i) / (ny - 1) * ylo +
(double) (i - 1) / (ny - 1) * yhi;
}
printf("\n");
printf("%s", "          X");
for (i = 1; i <= nx; ++i) {
printf("%8.2f", px[i - 1]);
}
printf("\n");
printf("%s", "     Y");
printf("\n");
for (i = ny; i >= 1; --i) {
for (j = 1; j <= nx; ++j) {
/* nag_2d_triang_eval (e01skc).
* A function to evaluate, at a set of points, the
* two-dimensional interpolant function generated by
* nag_2d_triang_interp (e01sjc).
*/
nag_2d_triang_eval(m, x, y, f, triang, grads, px[j - 1],
py[i - 1], &pf[j - 1], &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_2d_triang_eval (e01skc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
}
printf("%8.2f", py[i - 1]);
printf("%3s", "");
for (j = 1; j <= nx; ++j) {
printf("%8.2f", pf[j - 1]);
}
printf("\n");
}
END:
NAG_FREE(f);