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