1 /* interpolation/gsl_interp.h
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
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 /* Author: G. Jungman
21 */
22 #ifndef __GSL_INTERP_H__
23 #define __GSL_INTERP_H__
24 #include <stdlib.h>
25 #include "gsl_types.h"
26
27 #undef __BEGIN_DECLS
28 #undef __END_DECLS
29 #ifdef __cplusplus
30 # define __BEGIN_DECLS extern "C" {
31 # define __END_DECLS }
32 #else
33 # define __BEGIN_DECLS /* empty */
34 # define __END_DECLS /* empty */
35 #endif
36
37 __BEGIN_DECLS
38
39 /* evaluation accelerator */
40 typedef struct {
41 size_t cache; /* cache of index */
42 size_t miss_count; /* keep statistics */
43 size_t hit_count;
44 }
45 gsl_interp_accel;
46
47
48 /* interpolation object type */
49 typedef struct {
50 const char * name;
51 unsigned int min_size;
52 void * (*alloc) (size_t size);
53 int (*init) (void *, const double xa[], const double ya[], size_t size);
54 int (*eval) (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y);
55 int (*eval_deriv) (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y_p);
56 int (*eval_deriv2) (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y_pp);
57 int (*eval_integ) (const void *, const double xa[], const double ya[], size_t size, gsl_interp_accel *, double a, double b, double * result);
58 void (*free) (void *);
59
60 } gsl_interp_type;
61
62
63 /* general interpolation object */
64 typedef struct {
65 const gsl_interp_type * type;
66 double xmin;
67 double xmax;
68 size_t size;
69 void * state;
70 } gsl_interp;
71
72
73 /* available types */
74 GSL_VAR const gsl_interp_type * gsl_interp_linear;
75 GSL_VAR const gsl_interp_type * gsl_interp_polynomial;
76 GSL_VAR const gsl_interp_type * gsl_interp_cspline;
77 GSL_VAR const gsl_interp_type * gsl_interp_cspline_periodic;
78 GSL_VAR const gsl_interp_type * gsl_interp_akima;
79 GSL_VAR const gsl_interp_type * gsl_interp_akima_periodic;
80
81 gsl_interp_accel *
82 gsl_interp_accel_alloc(void);
83
84 size_t
85 gsl_interp_accel_find(gsl_interp_accel * a, const double x_array[], size_t size, double x);
86
87 int
88 gsl_interp_accel_reset (gsl_interp_accel * a);
89
90 void
91 gsl_interp_accel_free(gsl_interp_accel * a);
92
93 gsl_interp *
94 gsl_interp_alloc(const gsl_interp_type * T, size_t n);
95
96 int
97 gsl_interp_init(gsl_interp * obj, const double xa[], const double ya[], size_t size);
98
99 const char * gsl_interp_name(const gsl_interp * interp);
100 unsigned int gsl_interp_min_size(const gsl_interp * interp);
101
102
103 int
104 gsl_interp_eval_e(const gsl_interp * obj,
105 const double xa[], const double ya[], double x,
106 gsl_interp_accel * a, double * y);
107
108 double
109 gsl_interp_eval(const gsl_interp * obj,
110 const double xa[], const double ya[], double x,
111 gsl_interp_accel * a);
112
113 int
114 gsl_interp_eval_deriv_e(const gsl_interp * obj,
115 const double xa[], const double ya[], double x,
116 gsl_interp_accel * a,
117 double * d);
118
119 double
120 gsl_interp_eval_deriv(const gsl_interp * obj,
121 const double xa[], const double ya[], double x,
122 gsl_interp_accel * a);
123
124 int
125 gsl_interp_eval_deriv2_e(const gsl_interp * obj,
126 const double xa[], const double ya[], double x,
127 gsl_interp_accel * a,
128 double * d2);
129
130 double
131 gsl_interp_eval_deriv2(const gsl_interp * obj,
132 const double xa[], const double ya[], double x,
133 gsl_interp_accel * a);
134
135 int
136 gsl_interp_eval_integ_e(const gsl_interp * obj,
137 const double xa[], const double ya[],
138 double a, double b,
139 gsl_interp_accel * acc,
140 double * result);
141
142 double
143 gsl_interp_eval_integ(const gsl_interp * obj,
144 const double xa[], const double ya[],
145 double a, double b,
146 gsl_interp_accel * acc);
147
148 void
149 gsl_interp_free(gsl_interp * interp);
150
151 size_t gsl_interp_bsearch(const double x_array[], double x,
152 size_t index_lo, size_t index_hi);
153
154 #ifdef HAVE_INLINE
155 extern inline size_t
156 gsl_interp_bsearch(const double x_array[], double x,
157 size_t index_lo, size_t index_hi);
158
159 extern inline size_t
gsl_interp_bsearch(const double x_array[],double x,size_t index_lo,size_t index_hi)160 gsl_interp_bsearch(const double x_array[], double x,
161 size_t index_lo, size_t index_hi)
162 {
163 size_t ilo = index_lo;
164 size_t ihi = index_hi;
165 while(ihi > ilo + 1) {
166 size_t i = (ihi + ilo)/2;
167 if(x_array[i] > x)
168 ihi = i;
169 else
170 ilo = i;
171 }
172
173 return ilo;
174 }
175 #endif
176
177 #ifdef HAVE_INLINE
178 extern inline size_t
gsl_interp_accel_find(gsl_interp_accel * a,const double xa[],size_t len,double x)179 gsl_interp_accel_find(gsl_interp_accel * a, const double xa[], size_t len, double x)
180 {
181 size_t x_index = a->cache;
182
183 if(x < xa[x_index]) {
184 a->miss_count++;
185 a->cache = gsl_interp_bsearch(xa, x, 0, x_index);
186 }
187 else if(x > xa[x_index + 1]) {
188 a->miss_count++;
189 a->cache = gsl_interp_bsearch(xa, x, x_index, len-1);
190 }
191 else {
192 a->hit_count++;
193 }
194
195 return a->cache;
196 }
197 #endif /* HAVE_INLINE */
198
199
200 __END_DECLS
201
202 #endif /* __GSL_INTERP_H__ */
203