htswanalysis/MACS/lib/gsl/gsl-1.11/interpolation/linear.c

   1 /* interpolation/linear.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 #include <config.h>
  23 #include <stdlib.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_interp.h>
  26
  27 static int
  28 linear_init (void * vstate,
  29              const double x_array[],
  30              const double y_array[],
  31              size_t size)
  32 {
  33   return GSL_SUCCESS;
  34 }
  35
  36 static
  37 int
  38 linear_eval (const void * vstate,
  39              const double x_array[], const double y_array[], size_t size,
  40              double x,
  41              gsl_interp_accel * a,
  42              double *y)
  43 {
  44   double x_lo, x_hi;
  45   double y_lo, y_hi;
  46   double dx;
  47   size_t index;
  48
  49   if (a != 0)
  50     {
  51       index = gsl_interp_accel_find (a, x_array, size, x);
  52     }
  53   else
  54     {
  55       index = gsl_interp_bsearch (x_array, x, 0, size - 1);
  56     }
  57
  58   /* evaluate */
  59   x_lo = x_array[index];
  60   x_hi = x_array[index + 1];
  61   y_lo = y_array[index];
  62   y_hi = y_array[index + 1];
  63   dx = x_hi - x_lo;
  64   if (dx > 0.0)
  65     {
  66       *y = y_lo + (x - x_lo) / dx * (y_hi - y_lo);
  67       return GSL_SUCCESS;
  68     }
  69   else
  70     {
  71       *y = 0.0;
  72       return GSL_EINVAL;
  73     }
  74 }
  75
  76
  77 static
  78 int
  79 linear_eval_deriv (const void * vstate,
  80                    const double x_array[], const double y_array[], size_t size,
  81                    double x,
  82                    gsl_interp_accel * a,
  83                    double *dydx)
  84 {
  85   double x_lo, x_hi;
  86   double y_lo, y_hi;
  87   double dx;
  88   double dy;
  89   size_t index;
  90
  91   if (a != 0)
  92     {
  93       index = gsl_interp_accel_find (a, x_array, size, x);
  94     }
  95   else
  96     {
  97       index = gsl_interp_bsearch (x_array, x, 0, size - 1);
  98     }
  99
 100   /* evaluate */
 101   x_lo = x_array[index];
 102   x_hi = x_array[index + 1];
 103   y_lo = y_array[index];
 104   y_hi = y_array[index + 1];
 105   dx = x_hi - x_lo;
 106   dy = y_hi - y_lo;
 107   if (dx > 0.0)
 108     {
 109       *dydx = dy / dx;;
 110       return GSL_SUCCESS;
 111     }
 112   else
 113     {
 114       *dydx = 0.0;
 115       return GSL_EINVAL;
 116     }
 117 }
 118
 119
 120 static
 121 int
 122 linear_eval_deriv2 (const void * vstate,
 123                     const double x_array[], const double y_array[], size_t size,
 124                     double x,
 125                     gsl_interp_accel * a,
 126                     double *y_pp)
 127 {
 128   *y_pp = 0.0;
 129
 130   return GSL_SUCCESS;
 131 }
 132
 133
 134 static
 135 int
 136 linear_eval_integ (const void * vstate,
 137                    const double x_array[], const double y_array[], size_t size,
 138                    gsl_interp_accel * acc,
 139                    double a, double b,
 140                    double * result)
 141 {
 142   size_t i, index_a, index_b;
 143
 144   if (acc != 0)
 145     {
 146       index_a = gsl_interp_accel_find (acc, x_array, size, a);
 147       index_b = gsl_interp_accel_find (acc, x_array, size, b);
 148     }
 149   else
 150     {
 151       index_a = gsl_interp_bsearch (x_array, a, 0, size - 1);
 152       index_b = gsl_interp_bsearch (x_array, b, 0, size - 1);
 153     }
 154
 155     /* endpoints span more than one interval */
 156
 157   *result = 0.0;
 158
 159   /* interior intervals */
 160   for(i=index_a; i<=index_b; i++) {
 161     const double x_hi = x_array[i + 1];
 162     const double x_lo = x_array[i];
 163     const double y_lo = y_array[i];
 164     const double y_hi = y_array[i + 1];
 165     const double dx = x_hi - x_lo;
 166
 167     if(dx != 0.0) {
 168       if (i == index_a || i == index_b)
 169         {
 170           double x1 = (i == index_a) ? a : x_lo;
 171           double x2 = (i == index_b) ? b : x_hi;
 172           const double D = (y_hi-y_lo)/dx;
 173           *result += (x2-x1) * (y_lo + 0.5*D*((x2-x_lo)+(x1-x_lo)));
 174         }
 175       else
 176         {
 177           *result += 0.5 * dx * (y_lo + y_hi);
 178         }
 179     }
 180   }
 181
 182   return GSL_SUCCESS;
 183 }
 184
 185 static const gsl_interp_type linear_type =
 186 {
 187   "linear",
 188   2,
 189   NULL, /* alloc, not applicable */
 190   &linear_init,
 191   &linear_eval,
 192   &linear_eval_deriv,
 193   &linear_eval_deriv2,
 194   &linear_eval_integ,
 195   NULL, /* free, not applicable */
 196 };
 197
 198 const gsl_interp_type * gsl_interp_linear = &linear_type;