/* nag_2d_spline_interpolant (e01dac) Example Program. * * Copyright 1991 Numerical Algorithms Group. * * Mark 2, 1991. * * Mark 6 revised, 2000. * Mark 8 revised, 2004. */ #include #include #include #include #include #define F(I, J) f[my*(I)+(J)] #define FG(I, J) fg[npy*(I)+(J)] #define C(I, J) spline.c[my*(I)+(J)] int main(void) { Integer exit_status = 0, i, j, mx, my, npx, npy; NagError fail; Nag_2dSpline spline; double *f = 0, *fg = 0, step, *tx = 0, *ty = 0, *x = 0, xhi, xlo; double *y = 0, yhi, ylo; INIT_FAIL(fail); /* Initialise spline */ spline.lamda = 0; spline.mu = 0; spline.c = 0; printf( "nag_2d_spline_interpolant (e01dac) Example Program Results\n"); scanf("%*[^\n]"); /* Skip heading in data file */ /* Read the number of x points, mx, and the values of the * x co-ordinates. */ scanf("%ld%ld", &mx, &my); if (mx >= 4 && my >= 4) { if (!(f = NAG_ALLOC(mx*my, double)) || !(x = NAG_ALLOC(mx, double)) || !(y = NAG_ALLOC(my, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } } else { printf("Invalid mx or my.\n"); exit_status = 1; return exit_status; } for (i = 0; i < mx; i++) scanf("%lf", &x[i]); /* Read the number of y points, my, and the values of the * y co-ordinates. */ for (i = 0; i < my; i++) scanf("%lf", &y[i]); /* Read the function values at the grid points. */ for (j = 0; j < my; j++) for (i = 0; i < mx; i++) scanf("%lf", &F(i, j)); /* Generate the (x,y,f) interpolating bicubic B-spline. */ /* nag_2d_spline_interpolant (e01dac). * Interpolating function, bicubic spline interpolant, two * variables */ nag_2d_spline_interpolant(mx, my, x, y, f, &spline, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_2d_spline_interpolant (e01dac).\n%s\n", fail.message); exit_status = 1; goto END; } /* Print the knot sets, lamda and mu. */ printf("Distinct knots in x direction located at\n"); for (j = 3; j < spline.nx-3; j++) printf("%12.4f%s", spline.lamda[j], ((j-3)%5 == 4 || j == spline.nx-4)?"\n":" "); printf("\nDistinct knots in y direction located at\n"); for (j = 3; j < spline.ny-3; j++) printf("%12.4f%s", spline.mu[j], ((j-3)%5 == 4 || j == spline.ny-4)?"\n":" "); /* Print the spline coefficients. */ printf("\nThe B-Spline coefficients:\n"); for (i = 0; i < mx; i++) { for (j = 0; j < my; j++) printf("%9.4f", C(i, j)); printf("\n"); } /* Evaluate the spline on a regular rectangular grid at npx*npy * points over the domain (xlo to xhi) x (ylo to yhi). */ scanf("%ld%lf%lf", &npx, &xlo, &xhi); scanf("%ld%lf%lf", &npy, &ylo, &yhi); if (npx >= 1 && npy >= 1) { if (!(fg = NAG_ALLOC(npx*npy, double)) || !(tx = NAG_ALLOC(npx, double)) || !(ty = NAG_ALLOC(npy, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } } else { printf("Invalid npx or npy.\n"); exit_status = 1; return exit_status; } step = (xhi-xlo)/(double)(npx-1); printf("\nSpline evaluated on a regular mesh " " (x across, y down): \n "); /* Generate nx equispaced x co-ordinates. */ for (i = 0; i < npx; i++) { tx[i] = MIN(xlo+i*step, xhi); printf(" %5.2f ", tx[i]); } step = (yhi-ylo)/(npy-1); for (i = 0; i < npy; i++) ty[i] = MIN(ylo+i*step, yhi); /* Evaluate the spline. */ /* nag_2d_spline_eval_rect (e02dfc). * Evaluation of bicubic spline, at a mesh of points */ nag_2d_spline_eval_rect(npx, npy, tx, ty, fg, &spline, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_2d_spline_eval_rect (e02dfc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Print the results. */ printf("\n"); for (j = 0; j < npy; j++) { printf("%5.2f", ty[j]); for (i = 0; i < npx; i++) printf("%8.3f ", FG(i, j)); printf("\n"); } /* Free memory allocated by nag_2d_spline_interpolant (e01dac) */ END: NAG_FREE(spline.lamda); NAG_FREE(spline.mu); NAG_FREE(spline.c); if (f) NAG_FREE(f); if (x) NAG_FREE(x); if (y) NAG_FREE(y); if (fg) NAG_FREE(fg); if (tx) NAG_FREE(tx); if (ty) NAG_FREE(ty); return exit_status; }