1 /* interpolation/spline.c
2 *
3 * Copyright (C) 2001, 2007 Brian Gough
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 <string.h>
22 #include <gsl/gsl_errno.h>
23 #include <gsl/gsl_interp.h>
24 #include <gsl/gsl_spline.h>
25
26 gsl_spline *
gsl_spline_alloc(const gsl_interp_type * T,size_t size)27 gsl_spline_alloc (const gsl_interp_type * T, size_t size)
28 {
29 gsl_spline * spline = (gsl_spline *) malloc (sizeof(gsl_spline));
30
31 if (spline == NULL)
32 {
33 GSL_ERROR_NULL ("failed to allocate space for spline struct",
34 GSL_ENOMEM);
35 }
36
37 spline->interp = gsl_interp_alloc (T, size);
38
39 if (spline->interp == NULL)
40 {
41 free (spline);
42 GSL_ERROR_NULL ("failed to allocate space for interp", GSL_ENOMEM);
43 };
44
45 spline->x = (double *) malloc (size * sizeof(double));
46
47 if (spline->x == NULL)
48 {
49 gsl_interp_free(spline->interp);
50 free(spline);
51 GSL_ERROR_NULL ("failed to allocate space for x", GSL_ENOMEM);
52 }
53
54 spline->y = (double *) malloc (size * sizeof(double));
55
56 if (spline->y == NULL)
57 {
58 free(spline->x);
59 gsl_interp_free(spline->interp);
60 free(spline);
61 GSL_ERROR_NULL ("failed to allocate space for y", GSL_ENOMEM);
62 }
63
64 spline->size = size;
65
66 return spline;
67 }
68
69 int
gsl_spline_init(gsl_spline * spline,const double x_array[],const double y_array[],size_t size)70 gsl_spline_init (gsl_spline * spline, const double x_array[], const double y_array[], size_t size)
71 {
72 if (size != spline->size)
73 {
74 GSL_ERROR ("data must match size of spline object", GSL_EINVAL);
75 }
76
77 memcpy (spline->x, x_array, size * sizeof(double));
78 memcpy (spline->y, y_array, size * sizeof(double));
79
80 {
81 int status = gsl_interp_init (spline->interp, x_array, y_array, size);
82 return status;
83 }
84 }
85
86 const char *
gsl_spline_name(const gsl_spline * spline)87 gsl_spline_name(const gsl_spline * spline)
88 {
89 return gsl_interp_name(spline->interp);
90 }
91
92 unsigned int
gsl_spline_min_size(const gsl_spline * spline)93 gsl_spline_min_size(const gsl_spline * spline)
94 {
95 return gsl_interp_min_size(spline->interp);
96 }
97
98 void
gsl_spline_free(gsl_spline * spline)99 gsl_spline_free (gsl_spline * spline)
100 {
101 RETURN_IF_NULL (spline);
102 gsl_interp_free (spline->interp);
103 free (spline->x);
104 free (spline->y);
105 free (spline);
106 }
107
108 int
gsl_spline_eval_e(const gsl_spline * spline,double x,gsl_interp_accel * a,double * y)109 gsl_spline_eval_e (const gsl_spline * spline,
110 double x,
111 gsl_interp_accel * a, double *y)
112 {
113 return gsl_interp_eval_e (spline->interp,
114 spline->x, spline->y,
115 x, a, y);
116 }
117
118 double
gsl_spline_eval(const gsl_spline * spline,double x,gsl_interp_accel * a)119 gsl_spline_eval (const gsl_spline * spline,
120 double x,
121 gsl_interp_accel * a)
122 {
123 return gsl_interp_eval (spline->interp,
124 spline->x, spline->y,
125 x, a);
126 }
127
128
129 int
gsl_spline_eval_deriv_e(const gsl_spline * spline,double x,gsl_interp_accel * a,double * dydx)130 gsl_spline_eval_deriv_e (const gsl_spline * spline,
131 double x,
132 gsl_interp_accel * a,
133 double *dydx)
134 {
135 return gsl_interp_eval_deriv_e (spline->interp,
136 spline->x, spline->y,
137 x, a, dydx);
138 }
139
140 double
gsl_spline_eval_deriv(const gsl_spline * spline,double x,gsl_interp_accel * a)141 gsl_spline_eval_deriv (const gsl_spline * spline,
142 double x,
143 gsl_interp_accel * a)
144 {
145 return gsl_interp_eval_deriv (spline->interp,
146 spline->x, spline->y,
147 x, a);
148 }
149
150
151 int
gsl_spline_eval_deriv2_e(const gsl_spline * spline,double x,gsl_interp_accel * a,double * d2)152 gsl_spline_eval_deriv2_e (const gsl_spline * spline,
153 double x,
154 gsl_interp_accel * a,
155 double * d2)
156 {
157 return gsl_interp_eval_deriv2_e (spline->interp,
158 spline->x, spline->y,
159 x, a, d2);
160 }
161
162 double
gsl_spline_eval_deriv2(const gsl_spline * spline,double x,gsl_interp_accel * a)163 gsl_spline_eval_deriv2 (const gsl_spline * spline,
164 double x,
165 gsl_interp_accel * a)
166 {
167 return gsl_interp_eval_deriv2 (spline->interp,
168 spline->x, spline->y,
169 x, a);
170 }
171
172
173 int
gsl_spline_eval_integ_e(const gsl_spline * spline,double a,double b,gsl_interp_accel * acc,double * result)174 gsl_spline_eval_integ_e (const gsl_spline * spline,
175 double a, double b,
176 gsl_interp_accel * acc,
177 double * result)
178 {
179 return gsl_interp_eval_integ_e (spline->interp,
180 spline->x, spline->y,
181 a, b, acc, result);
182 }
183
184
185 double
gsl_spline_eval_integ(const gsl_spline * spline,double a,double b,gsl_interp_accel * acc)186 gsl_spline_eval_integ (const gsl_spline * spline,
187 double a, double b,
188 gsl_interp_accel * acc)
189 {
190 return gsl_interp_eval_integ (spline->interp,
191 spline->x, spline->y,
192 a, b, acc);
193 }
194
195