updated with main repo

1 files changed, 225 insertions, 0 deletions

diff --git a/thirdparty/includes/GSL/gsl/gsl_interp.h b/thirdparty/includes/GSL/gsl/gsl_interp.h
new file mode 100644
index 0000000..8321224
--- /dev/null
+++ b/thirdparty/includes/GSL/gsl/gsl_interp.h
@@ -0,0 +1,225 @@
+/* interpolation/gsl_interp.h
+ * 
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or (at
+ * your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+/* Author:  G. Jungman
+ */
+#ifndef __GSL_INTERP_H__
+#define __GSL_INTERP_H__
+#include <stdlib.h>
+#include <gsl/gsl_inline.h>
+#include <gsl/gsl_types.h>
+
+#undef __BEGIN_DECLS
+#undef __END_DECLS
+#ifdef __cplusplus
+# define __BEGIN_DECLS extern "C" {
+# define __END_DECLS }
+#else
+# define __BEGIN_DECLS /* empty */
+# define __END_DECLS /* empty */
+#endif
+
+__BEGIN_DECLS
+
+/* evaluation accelerator */
+typedef struct {
+  size_t  cache;        /* cache of index   */
+  size_t  miss_count;   /* keep statistics  */
+  size_t  hit_count;
+}
+gsl_interp_accel;
+
+
+/* interpolation object type */
+typedef struct {
+  const char * name;
+  unsigned int min_size;
+  void *  (*alloc) (size_t size);
+  int     (*init)    (void *, const double xa[], const double ya[], size_t size);
+  int     (*eval)    (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y);
+  int     (*eval_deriv)  (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y_p);
+  int     (*eval_deriv2) (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y_pp);
+  int     (*eval_integ)  (const void *, const double xa[], const double ya[], size_t size, gsl_interp_accel *, double a, double b, double * result);
+  void    (*free)         (void *);
+
+} gsl_interp_type;
+
+
+/* general interpolation object */
+typedef struct {
+  const gsl_interp_type * type;
+  double  xmin;
+  double  xmax;
+  size_t  size;
+  void * state;
+} gsl_interp;
+
+
+/* available types */
+GSL_VAR const gsl_interp_type * gsl_interp_linear;
+GSL_VAR const gsl_interp_type * gsl_interp_polynomial;
+GSL_VAR const gsl_interp_type * gsl_interp_cspline;
+GSL_VAR const gsl_interp_type * gsl_interp_cspline_periodic;
+GSL_VAR const gsl_interp_type * gsl_interp_akima;
+GSL_VAR const gsl_interp_type * gsl_interp_akima_periodic;
+GSL_VAR const gsl_interp_type * gsl_interp_steffen;
+
+gsl_interp_accel *
+gsl_interp_accel_alloc(void);
+
+int
+gsl_interp_accel_reset (gsl_interp_accel * a);
+
+void
+gsl_interp_accel_free(gsl_interp_accel * a);
+
+gsl_interp *
+gsl_interp_alloc(const gsl_interp_type * T, size_t n);
+     
+int
+gsl_interp_init(gsl_interp * obj, const double xa[], const double ya[], size_t size);
+
+const char * gsl_interp_name(const gsl_interp * interp);
+unsigned int gsl_interp_min_size(const gsl_interp * interp);
+unsigned int gsl_interp_type_min_size(const gsl_interp_type * T);
+
+
+int
+gsl_interp_eval_e(const gsl_interp * obj,
+                  const double xa[], const double ya[], double x,
+                  gsl_interp_accel * a, double * y);
+
+double
+gsl_interp_eval(const gsl_interp * obj,
+                const double xa[], const double ya[], double x,
+                gsl_interp_accel * a);
+
+int
+gsl_interp_eval_deriv_e(const gsl_interp * obj,
+                        const double xa[], const double ya[], double x,
+                        gsl_interp_accel * a,
+                        double * d);
+
+double
+gsl_interp_eval_deriv(const gsl_interp * obj,
+                      const double xa[], const double ya[], double x,
+                      gsl_interp_accel * a);
+
+int
+gsl_interp_eval_deriv2_e(const gsl_interp * obj,
+                         const double xa[], const double ya[], double x,
+                         gsl_interp_accel * a,
+                         double * d2);
+
+double
+gsl_interp_eval_deriv2(const gsl_interp * obj,
+                       const double xa[], const double ya[], double x,
+                       gsl_interp_accel * a);
+
+int
+gsl_interp_eval_integ_e(const gsl_interp * obj,
+                        const double xa[], const double ya[],
+                        double a, double b,
+                        gsl_interp_accel * acc,
+                        double * result);
+
+double
+gsl_interp_eval_integ(const gsl_interp * obj,
+                      const double xa[], const double ya[],
+                      double a, double b,
+                      gsl_interp_accel * acc);
+
+void
+gsl_interp_free(gsl_interp * interp);
+
+INLINE_DECL size_t
+gsl_interp_bsearch(const double x_array[], double x,
+                   size_t index_lo, size_t index_hi);
+
+#ifdef HAVE_INLINE
+
+/* Perform a binary search of an array of values.
+ * 
+ * The parameters index_lo and index_hi provide an initial bracket,
+ * and it is assumed that index_lo < index_hi. The resulting index
+ * is guaranteed to be strictly less than index_hi and greater than
+ * or equal to index_lo, so that the implicit bracket [index, index+1]
+ * always corresponds to a region within the implicit value range of
+ * the value array.
+ *
+ * Note that this means the relationship of 'x' to x_array[index]
+ * and x_array[index+1] depends on the result region, i.e. the
+ * behaviour at the boundaries may not correspond to what you
+ * expect. We have the following complete specification of the
+ * behaviour.
+ * Suppose the input is x_array[] = { x0, x1, ..., xN }
+ *    if ( x == x0 )           then  index == 0
+ *    if ( x > x0 && x <= x1 ) then  index == 0, and sim. for other interior pts
+ *    if ( x == xN )           then  index == N-1
+ *    if ( x > xN )            then  index == N-1
+ *    if ( x < x0 )            then  index == 0 
+ */
+
+INLINE_FUN size_t
+gsl_interp_bsearch(const double x_array[], double x,
+                   size_t index_lo, size_t index_hi)
+{
+  size_t ilo = index_lo;
+  size_t ihi = index_hi;
+  while(ihi > ilo + 1) {
+    size_t i = (ihi + ilo)/2;
+    if(x_array[i] > x)
+      ihi = i;
+    else
+      ilo = i;
+  }
+  
+  return ilo;
+}
+#endif
+
+INLINE_DECL size_t 
+gsl_interp_accel_find(gsl_interp_accel * a, const double x_array[], size_t size, double x);
+
+#ifdef HAVE_INLINE
+INLINE_FUN size_t
+gsl_interp_accel_find(gsl_interp_accel * a, const double xa[], size_t len, double x)
+{
+  size_t x_index = a->cache;
+ 
+  if(x < xa[x_index]) {
+    a->miss_count++;
+    a->cache = gsl_interp_bsearch(xa, x, 0, x_index);
+  }
+  else if(x >= xa[x_index + 1]) {
+    a->miss_count++;
+    a->cache = gsl_interp_bsearch(xa, x, x_index, len-1);
+  }
+  else {
+    a->hit_count++;
+  }
+  
+  return a->cache;
+}
+#endif /* HAVE_INLINE */
+
+
+__END_DECLS
+
+#endif /* __GSL_INTERP_H__ */