diff options
Diffstat (limited to 'gsl-1.9/interpolation/linear.c')
| -rw-r--r-- | gsl-1.9/interpolation/linear.c | 198 |
1 files changed, 198 insertions, 0 deletions
diff --git a/gsl-1.9/interpolation/linear.c b/gsl-1.9/interpolation/linear.c new file mode 100644 index 0000000..aa9ec00 --- /dev/null +++ b/gsl-1.9/interpolation/linear.c @@ -0,0 +1,198 @@ +/* interpolation/linear.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 2 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 + */ +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_interp.h> + +static int +linear_init (void * vstate, + const double x_array[], + const double y_array[], + size_t size) +{ + return GSL_SUCCESS; +} + +static +int +linear_eval (const void * vstate, + const double x_array[], const double y_array[], size_t size, + double x, + gsl_interp_accel * a, + double *y) +{ + double x_lo, x_hi; + double y_lo, y_hi; + double dx; + size_t index; + + if (a != 0) + { + index = gsl_interp_accel_find (a, x_array, size, x); + } + else + { + index = gsl_interp_bsearch (x_array, x, 0, size - 1); + } + + /* evaluate */ + x_lo = x_array[index]; + x_hi = x_array[index + 1]; + y_lo = y_array[index]; + y_hi = y_array[index + 1]; + dx = x_hi - x_lo; + if (dx > 0.0) + { + *y = y_lo + (x - x_lo) / dx * (y_hi - y_lo); + return GSL_SUCCESS; + } + else + { + *y = 0.0; + return GSL_EINVAL; + } +} + + +static +int +linear_eval_deriv (const void * vstate, + const double x_array[], const double y_array[], size_t size, + double x, + gsl_interp_accel * a, + double *dydx) +{ + double x_lo, x_hi; + double y_lo, y_hi; + double dx; + double dy; + size_t index; + + if (a != 0) + { + index = gsl_interp_accel_find (a, x_array, size, x); + } + else + { + index = gsl_interp_bsearch (x_array, x, 0, size - 1); + } + + /* evaluate */ + x_lo = x_array[index]; + x_hi = x_array[index + 1]; + y_lo = y_array[index]; + y_hi = y_array[index + 1]; + dx = x_hi - x_lo; + dy = y_hi - y_lo; + if (dx > 0.0) + { + *dydx = dy / dx;; + return GSL_SUCCESS; + } + else + { + *dydx = 0.0; + return GSL_EINVAL; + } +} + + +static +int +linear_eval_deriv2 (const void * vstate, + const double x_array[], const double y_array[], size_t size, + double x, + gsl_interp_accel * a, + double *y_pp) +{ + *y_pp = 0.0; + + return GSL_SUCCESS; +} + + +static +int +linear_eval_integ (const void * vstate, + const double x_array[], const double y_array[], size_t size, + gsl_interp_accel * acc, + double a, double b, + double * result) +{ + size_t i, index_a, index_b; + + if (acc != 0) + { + index_a = gsl_interp_accel_find (acc, x_array, size, a); + index_b = gsl_interp_accel_find (acc, x_array, size, b); + } + else + { + index_a = gsl_interp_bsearch (x_array, a, 0, size - 1); + index_b = gsl_interp_bsearch (x_array, b, 0, size - 1); + } + + /* endpoints span more than one interval */ + + *result = 0.0; + + /* interior intervals */ + for(i=index_a; i<=index_b; i++) { + const double x_hi = x_array[i + 1]; + const double x_lo = x_array[i]; + const double y_lo = y_array[i]; + const double y_hi = y_array[i + 1]; + const double dx = x_hi - x_lo; + + if(dx != 0.0) { + if (i == index_a || i == index_b) + { + double x1 = (i == index_a) ? a : x_lo; + double x2 = (i == index_b) ? b : x_hi; + const double D = (y_hi-y_lo)/dx; + *result += (x2-x1) * (y_lo + 0.5*D*((x2-x_lo)+(x1-x_lo))); + } + else + { + *result += 0.5 * dx * (y_lo + y_hi); + } + } + } + + return GSL_SUCCESS; +} + +static const gsl_interp_type linear_type = +{ + "linear", + 2, + NULL, /* alloc, not applicable */ + &linear_init, + &linear_eval, + &linear_eval_deriv, + &linear_eval_deriv2, + &linear_eval_integ, + NULL, /* free, not applicable */ +}; + +const gsl_interp_type * gsl_interp_linear = &linear_type; |