1 #include <stdio.h>
2 #include <stdlib.h>
3
4 #include <gsl/gsl_math.h>
5 #include <gsl/gsl_interp2d.h>
6 #include <gsl/gsl_spline2d.h>
7
8 int
main()9 main()
10 {
11 const gsl_interp2d_type *T = gsl_interp2d_bilinear;
12 const size_t N = 100; /* number of points to interpolate */
13 const double xa[] = { 0.0, 1.0 }; /* define unit square */
14 const double ya[] = { 0.0, 1.0 };
15 const size_t nx = sizeof(xa) / sizeof(double); /* x grid points */
16 const size_t ny = sizeof(ya) / sizeof(double); /* y grid points */
17 double *za = malloc(nx * ny * sizeof(double));
18 gsl_spline2d *spline = gsl_spline2d_alloc(T, nx, ny);
19 gsl_interp_accel *xacc = gsl_interp_accel_alloc();
20 gsl_interp_accel *yacc = gsl_interp_accel_alloc();
21 size_t i, j;
22
23 /* set z grid values */
24 gsl_spline2d_set(spline, za, 0, 0, 0.0);
25 gsl_spline2d_set(spline, za, 0, 1, 1.0);
26 gsl_spline2d_set(spline, za, 1, 1, 0.5);
27 gsl_spline2d_set(spline, za, 1, 0, 1.0);
28
29 /* initialize interpolation */
30 gsl_spline2d_init(spline, xa, ya, za, nx, ny);
31
32 /* interpolate N values in x and y and print out grid for plotting */
33 for (i = 0; i < N; ++i)
34 {
35 double xi = i / (N - 1.0);
36
37 for (j = 0; j < N; ++j)
38 {
39 double yj = j / (N - 1.0);
40 double zij = gsl_spline2d_eval(spline, xi, yj, xacc, yacc);
41
42 printf("%f %f %f\n", xi, yj, zij);
43 }
44 printf("\n");
45 }
46
47 gsl_spline2d_free(spline);
48 gsl_interp_accel_free(xacc);
49 gsl_interp_accel_free(yacc);
50 free(za);
51
52 return 0;
53 }
54