htswanalysis/MACS/lib/gsl/gsl-1.11/interpolation/gsl_interp.h

   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/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
 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
 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__ */