1 /* interpolation/interp.c
2 *
3 * Copyright (C) 2007 Brian Gough
4 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or (at
9 * your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 */
20
21 /* Author: G. Jungman
22 */
23 #include <config.h>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_math.h>
26 #include <gsl/gsl_interp.h>
27
28 #define DISCARD_STATUS(s) if ((s) != GSL_SUCCESS) { GSL_ERROR_VAL("interpolation error", (s), GSL_NAN); }
29
30 gsl_interp *
gsl_interp_alloc(const gsl_interp_type * T,size_t size)31 gsl_interp_alloc (const gsl_interp_type * T, size_t size)
32 {
33 gsl_interp * interp;
34
35 if (size < T->min_size)
36 {
37 GSL_ERROR_NULL ("insufficient number of points for interpolation type",
38 GSL_EINVAL);
39 }
40
41 interp = (gsl_interp *) malloc (sizeof(gsl_interp));
42
43 if (interp == NULL)
44 {
45 GSL_ERROR_NULL ("failed to allocate space for interp struct",
46 GSL_ENOMEM);
47 }
48
49 interp->type = T;
50 interp->size = size;
51
52 if (interp->type->alloc == NULL)
53 {
54 interp->state = NULL;
55 return interp;
56 }
57
58 interp->state = interp->type->alloc(size);
59
60 if (interp->state == NULL)
61 {
62 free (interp);
63 GSL_ERROR_NULL ("failed to allocate space for interp state", GSL_ENOMEM);
64 };
65
66 return interp;
67 }
68
69 int
gsl_interp_init(gsl_interp * interp,const double x_array[],const double y_array[],size_t size)70 gsl_interp_init (gsl_interp * interp, const double x_array[], const double y_array[], size_t size)
71 {
72 size_t i;
73
74 if (size != interp->size)
75 {
76 GSL_ERROR ("data must match size of interpolation object", GSL_EINVAL);
77 }
78
79 for (i = 1; i < size; i++)
80 {
81 if (!(x_array[i-1] < x_array[i]))
82 {
83 GSL_ERROR ("x values must be strictly increasing", GSL_EINVAL);
84 }
85 }
86
87 interp->xmin = x_array[0];
88 interp->xmax = x_array[size - 1];
89
90 {
91 int status = interp->type->init(interp->state, x_array, y_array, size);
92 return status;
93 }
94 }
95
96 const char *
gsl_interp_name(const gsl_interp * interp)97 gsl_interp_name(const gsl_interp * interp)
98 {
99 return interp->type->name;
100 }
101
102 unsigned int
gsl_interp_min_size(const gsl_interp * interp)103 gsl_interp_min_size(const gsl_interp * interp)
104 {
105 return interp->type->min_size;
106 }
107
108 unsigned int
gsl_interp_type_min_size(const gsl_interp_type * T)109 gsl_interp_type_min_size(const gsl_interp_type * T)
110 {
111 return T->min_size;
112 }
113
114 void
gsl_interp_free(gsl_interp * interp)115 gsl_interp_free (gsl_interp * interp)
116 {
117 RETURN_IF_NULL (interp);
118
119 if (interp->type->free)
120 interp->type->free (interp->state);
121 free (interp);
122 }
123
124
125
126 int
gsl_interp_eval_e(const gsl_interp * interp,const double xa[],const double ya[],double x,gsl_interp_accel * a,double * y)127 gsl_interp_eval_e (const gsl_interp * interp,
128 const double xa[], const double ya[], double x,
129 gsl_interp_accel * a, double *y)
130 {
131 if (x < interp->xmin || x > interp->xmax)
132 {
133 *y = GSL_NAN;
134 return GSL_EDOM;
135 }
136
137 return interp->type->eval (interp->state, xa, ya, interp->size, x, a, y);
138 }
139
140 double
gsl_interp_eval(const gsl_interp * interp,const double xa[],const double ya[],double x,gsl_interp_accel * a)141 gsl_interp_eval (const gsl_interp * interp,
142 const double xa[], const double ya[], double x,
143 gsl_interp_accel * a)
144 {
145 double y;
146 int status;
147
148 if (x < interp->xmin || x > interp->xmax)
149 {
150 GSL_ERROR_VAL("interpolation error", GSL_EDOM, GSL_NAN);
151 }
152
153 status = interp->type->eval (interp->state, xa, ya, interp->size, x, a, &y);
154
155 DISCARD_STATUS(status);
156
157 return y;
158 }
159
160
161 int
gsl_interp_eval_deriv_e(const gsl_interp * interp,const double xa[],const double ya[],double x,gsl_interp_accel * a,double * dydx)162 gsl_interp_eval_deriv_e (const gsl_interp * interp,
163 const double xa[], const double ya[], double x,
164 gsl_interp_accel * a,
165 double *dydx)
166 {
167 if (x < interp->xmin || x > interp->xmax)
168 {
169 *dydx = GSL_NAN;
170 return GSL_EDOM;
171 }
172
173 return interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, dydx);
174 }
175
176 double
gsl_interp_eval_deriv(const gsl_interp * interp,const double xa[],const double ya[],double x,gsl_interp_accel * a)177 gsl_interp_eval_deriv (const gsl_interp * interp,
178 const double xa[], const double ya[], double x,
179 gsl_interp_accel * a)
180 {
181 double dydx;
182 int status;
183
184 if (x < interp->xmin || x > interp->xmax)
185 {
186 GSL_ERROR_VAL("interpolation error", GSL_EDOM, GSL_NAN);
187 }
188
189 status = interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, &dydx);
190
191 DISCARD_STATUS(status);
192
193 return dydx;
194 }
195
196
197 int
gsl_interp_eval_deriv2_e(const gsl_interp * interp,const double xa[],const double ya[],double x,gsl_interp_accel * a,double * d2)198 gsl_interp_eval_deriv2_e (const gsl_interp * interp,
199 const double xa[], const double ya[], double x,
200 gsl_interp_accel * a,
201 double * d2)
202 {
203 if (x < interp->xmin || x > interp->xmax)
204 {
205 *d2 = GSL_NAN;
206 return GSL_EDOM;
207 }
208
209 return interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, d2);
210 }
211
212 double
gsl_interp_eval_deriv2(const gsl_interp * interp,const double xa[],const double ya[],double x,gsl_interp_accel * a)213 gsl_interp_eval_deriv2 (const gsl_interp * interp,
214 const double xa[], const double ya[], double x,
215 gsl_interp_accel * a)
216 {
217 double d2;
218 int status;
219
220 if (x < interp->xmin || x > interp->xmax)
221 {
222 GSL_ERROR_VAL("interpolation error", GSL_EDOM, GSL_NAN);
223 }
224
225 status = interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, &d2);
226
227 DISCARD_STATUS(status);
228
229 return d2;
230 }
231
232
233 int
gsl_interp_eval_integ_e(const gsl_interp * interp,const double xa[],const double ya[],double a,double b,gsl_interp_accel * acc,double * result)234 gsl_interp_eval_integ_e (const gsl_interp * interp,
235 const double xa[], const double ya[],
236 double a, double b,
237 gsl_interp_accel * acc,
238 double * result)
239 {
240 if (a > b || a < interp->xmin || b > interp->xmax)
241 {
242 *result = GSL_NAN;
243 return GSL_EDOM;
244 }
245 else if(a == b)
246 {
247 *result = 0.0;
248 return GSL_SUCCESS;
249 }
250
251 return interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, result);
252 }
253
254
255 double
gsl_interp_eval_integ(const gsl_interp * interp,const double xa[],const double ya[],double a,double b,gsl_interp_accel * acc)256 gsl_interp_eval_integ (const gsl_interp * interp,
257 const double xa[], const double ya[],
258 double a, double b,
259 gsl_interp_accel * acc)
260 {
261 double result;
262 int status;
263
264 if (a > b || a < interp->xmin || b > interp->xmax)
265 {
266 GSL_ERROR_VAL("interpolation error", GSL_EDOM, GSL_NAN);
267 }
268 else if(a == b)
269 {
270 return 0.0;
271 }
272
273 status = interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, &result);
274
275 DISCARD_STATUS(status);
276
277 return result;
278 }
279
280
281