include/gsl/gsl_interp.h

0001 /* interpolation/gsl_interp.h
0002  *
0003  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
0004  *
0005  * This program is free software; you can redistribute it and/or modify
0006  * it under the terms of the GNU General Public License as published by
0007  * the Free Software Foundation; either version 3 of the License, or (at
0008  * your option) any later version.
0009  *
0010  * This program is distributed in the hope that it will be useful, but
0011  * WITHOUT ANY WARRANTY; without even the implied warranty of
0012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0013  * General Public License for more details.
0014  *
0015  * You should have received a copy of the GNU General Public License
0016  * along with this program; if not, write to the Free Software
0017  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
0018  */
0019
0020 /* Author:  G. Jungman
0021  */
0022 #ifndef __GSL_INTERP_H__
0023 #define __GSL_INTERP_H__
0024 #include <stdlib.h>
0025 #include <gsl/gsl_inline.h>
0026 #include <gsl/gsl_types.h>
0027
0028 #undef __BEGIN_DECLS
0029 #undef __END_DECLS
0030 #ifdef __cplusplus
0031 # define __BEGIN_DECLS extern "C" {
0032 # define __END_DECLS }
0033 #else
0034 # define __BEGIN_DECLS /* empty */
0035 # define __END_DECLS /* empty */
0036 #endif
0037
0038 __BEGIN_DECLS
0039
0040 /* evaluation accelerator */
0041 typedef struct {
0042   size_t  cache;        /* cache of index   */
0043   size_t  miss_count;   /* keep statistics  */
0044   size_t  hit_count;
0045 }
0046 gsl_interp_accel;
0047
0048
0049 /* interpolation object type */
0050 typedef struct {
0051   const char * name;
0052   unsigned int min_size;
0053   void *  (*alloc) (size_t size);
0054   int     (*init)    (void *, const double xa[], const double ya[], size_t size);
0055   int     (*eval)    (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y);
0056   int     (*eval_deriv)  (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y_p);
0057   int     (*eval_deriv2) (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y_pp);
0058   int     (*eval_integ)  (const void *, const double xa[], const double ya[], size_t size, gsl_interp_accel *, double a, double b, double * result);
0059   void    (*free)         (void *);
0060
0061 } gsl_interp_type;
0062
0063
0064 /* general interpolation object */
0065 typedef struct {
0066   const gsl_interp_type * type;
0067   double  xmin;
0068   double  xmax;
0069   size_t  size;
0070   void * state;
0071 } gsl_interp;
0072
0073
0074 /* available types */
0075 GSL_VAR const gsl_interp_type * gsl_interp_linear;
0076 GSL_VAR const gsl_interp_type * gsl_interp_polynomial;
0077 GSL_VAR const gsl_interp_type * gsl_interp_cspline;
0078 GSL_VAR const gsl_interp_type * gsl_interp_cspline_periodic;
0079 GSL_VAR const gsl_interp_type * gsl_interp_akima;
0080 GSL_VAR const gsl_interp_type * gsl_interp_akima_periodic;
0081 GSL_VAR const gsl_interp_type * gsl_interp_steffen;
0082
0083 gsl_interp_accel *
0084 gsl_interp_accel_alloc(void);
0085
0086 int
0087 gsl_interp_accel_reset (gsl_interp_accel * a);
0088
0089 void
0090 gsl_interp_accel_free(gsl_interp_accel * a);
0091
0092 gsl_interp *
0093 gsl_interp_alloc(const gsl_interp_type * T, size_t n);
0094
0095 int
0096 gsl_interp_init(gsl_interp * obj, const double xa[], const double ya[], size_t size);
0097
0098 const char * gsl_interp_name(const gsl_interp * interp);
0099 unsigned int gsl_interp_min_size(const gsl_interp * interp);
0100 unsigned int gsl_interp_type_min_size(const gsl_interp_type * T);
0101
0102
0103 int
0104 gsl_interp_eval_e(const gsl_interp * obj,
0105                   const double xa[], const double ya[], double x,
0106                   gsl_interp_accel * a, double * y);
0107
0108 double
0109 gsl_interp_eval(const gsl_interp * obj,
0110                 const double xa[], const double ya[], double x,
0111                 gsl_interp_accel * a);
0112
0113 int
0114 gsl_interp_eval_deriv_e(const gsl_interp * obj,
0115                         const double xa[], const double ya[], double x,
0116                         gsl_interp_accel * a,
0117                         double * d);
0118
0119 double
0120 gsl_interp_eval_deriv(const gsl_interp * obj,
0121                       const double xa[], const double ya[], double x,
0122                       gsl_interp_accel * a);
0123
0124 int
0125 gsl_interp_eval_deriv2_e(const gsl_interp * obj,
0126                          const double xa[], const double ya[], double x,
0127                          gsl_interp_accel * a,
0128                          double * d2);
0129
0130 double
0131 gsl_interp_eval_deriv2(const gsl_interp * obj,
0132                        const double xa[], const double ya[], double x,
0133                        gsl_interp_accel * a);
0134
0135 int
0136 gsl_interp_eval_integ_e(const gsl_interp * obj,
0137                         const double xa[], const double ya[],
0138                         double a, double b,
0139                         gsl_interp_accel * acc,
0140                         double * result);
0141
0142 double
0143 gsl_interp_eval_integ(const gsl_interp * obj,
0144                       const double xa[], const double ya[],
0145                       double a, double b,
0146                       gsl_interp_accel * acc);
0147
0148 void
0149 gsl_interp_free(gsl_interp * interp);
0150
0151 INLINE_DECL size_t
0152 gsl_interp_bsearch(const double x_array[], double x,
0153                    size_t index_lo, size_t index_hi);
0154
0155 #ifdef HAVE_INLINE
0156
0157 /* Perform a binary search of an array of values.
0158  *
0159  * The parameters index_lo and index_hi provide an initial bracket,
0160  * and it is assumed that index_lo < index_hi. The resulting index
0161  * is guaranteed to be strictly less than index_hi and greater than
0162  * or equal to index_lo, so that the implicit bracket [index, index+1]
0163  * always corresponds to a region within the implicit value range of
0164  * the value array.
0165  *
0166  * Note that this means the relationship of 'x' to x_array[index]
0167  * and x_array[index+1] depends on the result region, i.e. the
0168  * behaviour at the boundaries may not correspond to what you
0169  * expect. We have the following complete specification of the
0170  * behaviour.
0171  * Suppose the input is x_array[] = { x0, x1, ..., xN }
0172  *    if ( x == x0 )           then  index == 0
0173  *    if ( x > x0 && x <= x1 ) then  index == 0, and sim. for other interior pts
0174  *    if ( x == xN )           then  index == N-1
0175  *    if ( x > xN )            then  index == N-1
0176  *    if ( x < x0 )            then  index == 0
0177  */
0178
0179 INLINE_FUN size_t
0180 gsl_interp_bsearch(const double x_array[], double x,
0181                    size_t index_lo, size_t index_hi)
0182 {
0183   size_t ilo = index_lo;
0184   size_t ihi = index_hi;
0185   while(ihi > ilo + 1) {
0186     size_t i = (ihi + ilo)/2;
0187     if(x_array[i] > x)
0188       ihi = i;
0189     else
0190       ilo = i;
0191   }
0192
0193   return ilo;
0194 }
0195 #endif
0196
0197 INLINE_DECL size_t
0198 gsl_interp_accel_find(gsl_interp_accel * a, const double x_array[], size_t size, double x);
0199
0200 #ifdef HAVE_INLINE
0201 INLINE_FUN size_t
0202 gsl_interp_accel_find(gsl_interp_accel * a, const double xa[], size_t len, double x)
0203 {
0204   size_t x_index = a->cache;
0205
0206   if(x < xa[x_index]) {
0207     a->miss_count++;
0208     a->cache = gsl_interp_bsearch(xa, x, 0, x_index);
0209   }
0210   else if(x >= xa[x_index + 1]) {
0211     a->miss_count++;
0212     a->cache = gsl_interp_bsearch(xa, x, x_index, len-1);
0213   }
0214   else {
0215     a->hit_count++;
0216   }
0217
0218   return a->cache;
0219 }
0220 #endif /* HAVE_INLINE */
0221
0222
0223 __END_DECLS
0224
0225 #endif /* __GSL_INTERP_H__ */