1 /* interpolation/interp2d.c
2 *
3 * Copyright 2012 David Zaslavsky
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include <config.h>
21 #include <stdlib.h>
22 #include <gsl/gsl_errno.h>
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_interp.h>
25 #include <gsl/gsl_interp2d.h>
26
27 /**
28 * Triggers a GSL error if the argument is not equal to GSL_SUCCESS.
29 * If the argument is GSL_SUCCESS, this does nothing.
30 */
31 #define DISCARD_STATUS(s) if ((s) != GSL_SUCCESS) { GSL_ERROR_VAL("interpolation error", (s), GSL_NAN); }
32
33 #define IDX2D(i, j, w) ((j) * ((w)->xsize) + (i))
34
35 gsl_interp2d *
gsl_interp2d_alloc(const gsl_interp2d_type * T,const size_t xsize,const size_t ysize)36 gsl_interp2d_alloc(const gsl_interp2d_type * T, const size_t xsize,
37 const size_t ysize)
38 {
39 gsl_interp2d * interp;
40
41 if (xsize < T->min_size || ysize < T->min_size)
42 {
43 GSL_ERROR_NULL ("insufficient number of points for interpolation type",
44 GSL_EINVAL);
45 }
46
47 interp = (gsl_interp2d *) calloc(1, sizeof(gsl_interp2d));
48 if (interp == NULL)
49 {
50 GSL_ERROR_NULL ("failed to allocate space for gsl_interp2d struct",
51 GSL_ENOMEM);
52 }
53
54 interp->type = T;
55 interp->xsize = xsize;
56 interp->ysize = ysize;
57
58 if (interp->type->alloc == NULL)
59 {
60 interp->state = NULL;
61 return interp;
62 }
63
64 interp->state = interp->type->alloc(xsize, ysize);
65 if (interp->state == NULL)
66 {
67 free(interp);
68 GSL_ERROR_NULL ("failed to allocate space for gsl_interp2d state",
69 GSL_ENOMEM);
70 }
71
72 return interp;
73 } /* gsl_interp2d_alloc() */
74
75 void
gsl_interp2d_free(gsl_interp2d * interp)76 gsl_interp2d_free (gsl_interp2d * interp)
77 {
78 RETURN_IF_NULL(interp);
79
80 if (interp->type->free)
81 interp->type->free(interp->state);
82
83 free(interp);
84 } /* gsl_interp2d_free() */
85
86 int
gsl_interp2d_init(gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const size_t xsize,const size_t ysize)87 gsl_interp2d_init (gsl_interp2d * interp, const double xarr[], const double yarr[],
88 const double zarr[], const size_t xsize, const size_t ysize)
89 {
90 size_t i;
91
92 if (xsize != interp->xsize || ysize != interp->ysize)
93 {
94 GSL_ERROR("data must match size of interpolation object", GSL_EINVAL);
95 }
96
97 for (i = 1; i < xsize; i++)
98 {
99 if (xarr[i-1] >= xarr[i])
100 {
101 GSL_ERROR("x values must be strictly increasing", GSL_EINVAL);
102 }
103 }
104
105 for (i = 1; i < ysize; i++)
106 {
107 if (yarr[i-1] >= yarr[i])
108 {
109 GSL_ERROR("y values must be strictly increasing", GSL_EINVAL);
110 }
111 }
112
113 interp->xmin = xarr[0];
114 interp->xmax = xarr[xsize - 1];
115 interp->ymin = yarr[0];
116 interp->ymax = yarr[ysize - 1];
117
118 {
119 int status = interp->type->init(interp->state, xarr, yarr, zarr,
120 xsize, ysize);
121 return status;
122 }
123 } /* gsl_interp2d_init() */
124
125 /*
126 * A wrapper function that checks boundary conditions, calls an evaluator
127 * which implements the actual calculation of the function value or
128 * derivative etc., and checks the return status.
129 */
130 static int
interp2d_eval(int (* evaluator)(const void *,const double xa[],const double ya[],const double za[],size_t xsize,size_t ysize,double x,double y,gsl_interp_accel *,gsl_interp_accel *,double * z),const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya,double * result)131 interp2d_eval(int (*evaluator)(const void *, const double xa[], const double ya[],
132 const double za[], size_t xsize, size_t ysize,
133 double x, double y, gsl_interp_accel *,
134 gsl_interp_accel *, double * z),
135 const gsl_interp2d * interp, const double xarr[],
136 const double yarr[], const double zarr[],
137 const double x, const double y,
138 gsl_interp_accel * xa, gsl_interp_accel * ya,
139 double * result)
140 {
141 if (x < interp->xmin || x > interp->xmax)
142 {
143 GSL_ERROR ("interpolation x value out of range", GSL_EDOM);
144 }
145 else if (y < interp->ymin || y > interp->ymax)
146 {
147 GSL_ERROR ("interpolation y value out of range", GSL_EDOM);
148 }
149
150 return evaluator(interp->state, xarr, yarr, zarr,
151 interp->xsize, interp->ysize,
152 x, y, xa, ya, result);
153 }
154
155 /*
156 * Another wrapper function that serves as a drop-in replacement for
157 * interp2d_eval but does not check the bounds. This can be used
158 * for extrapolation.
159 */
160 static int
interp2d_eval_extrap(int (* evaluator)(const void *,const double xa[],const double ya[],const double za[],size_t xsize,size_t ysize,double x,double y,gsl_interp_accel *,gsl_interp_accel *,double * z),const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya,double * result)161 interp2d_eval_extrap(int (*evaluator)(const void *, const double xa[],
162 const double ya[], const double za[],
163 size_t xsize, size_t ysize,
164 double x, double y,
165 gsl_interp_accel *,
166 gsl_interp_accel *, double * z),
167 const gsl_interp2d * interp, const double xarr[],
168 const double yarr[], const double zarr[],
169 const double x, const double y,
170 gsl_interp_accel * xa, gsl_interp_accel * ya,
171 double * result)
172 {
173 return evaluator(interp->state, xarr, yarr, zarr,
174 interp->xsize, interp->ysize, x, y, xa, ya, result);
175 }
176
177 double
gsl_interp2d_eval(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya)178 gsl_interp2d_eval (const gsl_interp2d * interp, const double xarr[],
179 const double yarr[], const double zarr[],
180 const double x, const double y,
181 gsl_interp_accel * xa, gsl_interp_accel * ya)
182 {
183 double z;
184 int status = gsl_interp2d_eval_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
185 DISCARD_STATUS(status)
186 return z;
187 } /* gsl_interp2d_eval() */
188
189 double
gsl_interp2d_eval_extrap(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya)190 gsl_interp2d_eval_extrap (const gsl_interp2d * interp,
191 const double xarr[],
192 const double yarr[],
193 const double zarr[],
194 const double x,
195 const double y,
196 gsl_interp_accel * xa,
197 gsl_interp_accel * ya)
198 {
199 double z;
200 int status =
201 interp2d_eval_extrap(interp->type->eval, interp,
202 xarr, yarr, zarr, x, y, xa, ya, &z);
203 DISCARD_STATUS(status)
204 return z;
205 }
206
207 int
gsl_interp2d_eval_e(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya,double * z)208 gsl_interp2d_eval_e (const gsl_interp2d * interp, const double xarr[],
209 const double yarr[], const double zarr[],
210 const double x, const double y,
211 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
212 {
213 return interp2d_eval(interp->type->eval, interp,
214 xarr, yarr, zarr, x, y, xa, ya, z);
215 } /* gsl_interp2d_eval_e() */
216
217 #ifndef GSL_DISABLE_DEPRECATED
218
219 int
gsl_interp2d_eval_e_extrap(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya,double * z)220 gsl_interp2d_eval_e_extrap (const gsl_interp2d * interp,
221 const double xarr[], const double yarr[],
222 const double zarr[], const double x,
223 const double y, gsl_interp_accel * xa,
224 gsl_interp_accel * ya, double * z)
225 {
226 return interp2d_eval_extrap(interp->type->eval, interp,
227 xarr, yarr, zarr, x, y, xa, ya, z);
228 }
229
230 #endif /* !GSL_DISABLE_DEPRECATED */
231
232 int
gsl_interp2d_eval_extrap_e(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya,double * z)233 gsl_interp2d_eval_extrap_e (const gsl_interp2d * interp,
234 const double xarr[], const double yarr[],
235 const double zarr[], const double x,
236 const double y, gsl_interp_accel * xa,
237 gsl_interp_accel * ya, double * z)
238 {
239 return interp2d_eval_extrap(interp->type->eval, interp,
240 xarr, yarr, zarr, x, y, xa, ya, z);
241 }
242
243 double
gsl_interp2d_eval_deriv_x(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya)244 gsl_interp2d_eval_deriv_x (const gsl_interp2d * interp, const double xarr[],
245 const double yarr[], const double zarr[],
246 const double x, const double y,
247 gsl_interp_accel * xa, gsl_interp_accel * ya)
248 {
249 double z;
250 int status = gsl_interp2d_eval_deriv_x_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
251 DISCARD_STATUS(status)
252 return z;
253 }
254
255 int
gsl_interp2d_eval_deriv_x_e(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya,double * z)256 gsl_interp2d_eval_deriv_x_e (const gsl_interp2d * interp, const double xarr[],
257 const double yarr[], const double zarr[],
258 const double x, const double y,
259 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
260 {
261 return interp2d_eval(interp->type->eval_deriv_x, interp,
262 xarr, yarr, zarr, x, y, xa, ya, z);
263 }
264
265 double
gsl_interp2d_eval_deriv_y(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya)266 gsl_interp2d_eval_deriv_y (const gsl_interp2d * interp, const double xarr[],
267 const double yarr[], const double zarr[],
268 const double x, const double y,
269 gsl_interp_accel * xa, gsl_interp_accel * ya)
270 {
271 double z;
272 int status = gsl_interp2d_eval_deriv_y_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
273 DISCARD_STATUS(status)
274 return z;
275 }
276
277 int
gsl_interp2d_eval_deriv_y_e(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya,double * z)278 gsl_interp2d_eval_deriv_y_e (const gsl_interp2d * interp, const double xarr[],
279 const double yarr[], const double zarr[],
280 const double x, const double y,
281 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
282 {
283 return interp2d_eval(interp->type->eval_deriv_y, interp,
284 xarr, yarr, zarr, x, y, xa, ya, z);
285 }
286
287 double
gsl_interp2d_eval_deriv_xx(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya)288 gsl_interp2d_eval_deriv_xx (const gsl_interp2d * interp, const double xarr[],
289 const double yarr[], const double zarr[],
290 const double x, const double y,
291 gsl_interp_accel * xa, gsl_interp_accel * ya)
292 {
293 double z;
294 int status = gsl_interp2d_eval_deriv_xx_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
295 DISCARD_STATUS(status)
296 return z;
297 }
298
299 int
gsl_interp2d_eval_deriv_xx_e(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya,double * z)300 gsl_interp2d_eval_deriv_xx_e (const gsl_interp2d * interp, const double xarr[],
301 const double yarr[], const double zarr[],
302 const double x, const double y,
303 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
304 {
305 return interp2d_eval(interp->type->eval_deriv_xx, interp,
306 xarr, yarr, zarr, x, y, xa, ya, z);
307 }
308
309 double
gsl_interp2d_eval_deriv_yy(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya)310 gsl_interp2d_eval_deriv_yy (const gsl_interp2d * interp, const double xarr[],
311 const double yarr[], const double zarr[],
312 const double x, const double y,
313 gsl_interp_accel * xa, gsl_interp_accel * ya)
314 {
315 double z;
316 int status = gsl_interp2d_eval_deriv_yy_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
317 DISCARD_STATUS(status)
318 return z;
319 }
320
321 int
gsl_interp2d_eval_deriv_yy_e(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya,double * z)322 gsl_interp2d_eval_deriv_yy_e (const gsl_interp2d * interp, const double xarr[],
323 const double yarr[], const double zarr[],
324 const double x, const double y,
325 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
326 {
327 return interp2d_eval(interp->type->eval_deriv_yy, interp,
328 xarr, yarr, zarr, x, y, xa, ya, z);
329 }
330
331 double
gsl_interp2d_eval_deriv_xy(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya)332 gsl_interp2d_eval_deriv_xy (const gsl_interp2d * interp, const double xarr[],
333 const double yarr[], const double zarr[],
334 const double x, const double y,
335 gsl_interp_accel * xa, gsl_interp_accel * ya)
336 {
337 double z;
338 int status = gsl_interp2d_eval_deriv_xy_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
339 DISCARD_STATUS(status)
340 return z;
341 }
342
343 int
gsl_interp2d_eval_deriv_xy_e(const gsl_interp2d * interp,const double xarr[],const double yarr[],const double zarr[],const double x,const double y,gsl_interp_accel * xa,gsl_interp_accel * ya,double * z)344 gsl_interp2d_eval_deriv_xy_e (const gsl_interp2d * interp, const double xarr[],
345 const double yarr[], const double zarr[],
346 const double x, const double y,
347 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
348 {
349 return interp2d_eval(interp->type->eval_deriv_xy, interp,
350 xarr, yarr, zarr, x, y, xa, ya, z);
351 }
352
353 size_t
gsl_interp2d_type_min_size(const gsl_interp2d_type * T)354 gsl_interp2d_type_min_size(const gsl_interp2d_type * T)
355 {
356 return T->min_size;
357 }
358
359 size_t
gsl_interp2d_min_size(const gsl_interp2d * interp)360 gsl_interp2d_min_size(const gsl_interp2d * interp)
361 {
362 return interp->type->min_size;
363 }
364
365 const char *
gsl_interp2d_name(const gsl_interp2d * interp)366 gsl_interp2d_name(const gsl_interp2d * interp)
367 {
368 return interp->type->name;
369 }
370
371 size_t
gsl_interp2d_idx(const gsl_interp2d * interp,const size_t i,const size_t j)372 gsl_interp2d_idx(const gsl_interp2d * interp,
373 const size_t i, const size_t j)
374 {
375 if (i >= interp->xsize)
376 {
377 GSL_ERROR_VAL ("x index out of range", GSL_ERANGE, 0);
378 }
379 else if (j >= interp->ysize)
380 {
381 GSL_ERROR_VAL ("y index out of range", GSL_ERANGE, 0);
382 }
383 else
384 {
385 return IDX2D(i, j, interp);
386 }
387 } /* gsl_interp2d_idx() */
388
389 int
gsl_interp2d_set(const gsl_interp2d * interp,double zarr[],const size_t i,const size_t j,const double z)390 gsl_interp2d_set(const gsl_interp2d * interp, double zarr[],
391 const size_t i, const size_t j, const double z)
392 {
393 if (i >= interp->xsize)
394 {
395 GSL_ERROR ("x index out of range", GSL_ERANGE);
396 }
397 else if (j >= interp->ysize)
398 {
399 GSL_ERROR ("y index out of range", GSL_ERANGE);
400 }
401 else
402 {
403 zarr[IDX2D(i, j, interp)] = z;
404 return GSL_SUCCESS;
405 }
406 } /* gsl_interp2d_set() */
407
408 double
gsl_interp2d_get(const gsl_interp2d * interp,const double zarr[],const size_t i,const size_t j)409 gsl_interp2d_get(const gsl_interp2d * interp, const double zarr[],
410 const size_t i, const size_t j)
411 {
412 if (i >= interp->xsize)
413 {
414 GSL_ERROR_VAL ("x index out of range", GSL_ERANGE, 0);
415 }
416 else if (j >= interp->ysize)
417 {
418 GSL_ERROR_VAL ("y index out of range", GSL_ERANGE, 0);
419 }
420 else
421 {
422 return zarr[IDX2D(i, j, interp)];
423 }
424 } /* gsl_interp2d_get() */
425
426 #undef IDX2D
427