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

   1 /* interpolation/interp.c
   2  *
   3  * Copyright (C) 2007 Brian Gough
   4  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19  */
  20
  21 /* Author:  G. Jungman
  22  */
  23 #include <config.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_math.h>
  26 #include <gsl/gsl_interp.h>
  27
  28 #define DISCARD_STATUS(s) if ((s) != GSL_SUCCESS) { GSL_ERROR_VAL("interpolation error", (s),  GSL_NAN); }
  29
  30 gsl_interp *
  31 gsl_interp_alloc (const gsl_interp_type * T, size_t size)
  32 {
  33   gsl_interp * interp;
  34
  35   if (size < T->min_size)
  36     {
  37       GSL_ERROR_NULL ("insufficient number of points for interpolation type",
  38                       GSL_EINVAL);
  39     }
  40
  41   interp = (gsl_interp *) malloc (sizeof(gsl_interp));
  42
  43   if (interp == NULL)
  44     {
  45       GSL_ERROR_NULL ("failed to allocate space for interp struct",
  46                       GSL_ENOMEM);
  47     }
  48
  49   interp->type = T;
  50   interp->size = size;
  51
  52   if (interp->type->alloc == NULL)
  53     {
  54       interp->state = NULL;
  55       return interp;
  56     }
  57
  58   interp->state = interp->type->alloc(size);
  59
  60   if (interp->state == NULL)
  61     {
  62       free (interp);
  63       GSL_ERROR_NULL ("failed to allocate space for interp state", GSL_ENOMEM);
  64     };
  65
  66   return interp;
  67 }
  68
  69 int
  70 gsl_interp_init (gsl_interp * interp, const double x_array[], const double y_array[], size_t size)
  71 {
  72   size_t i;
  73
  74   if (size != interp->size)
  75     {
  76       GSL_ERROR ("data must match size of interpolation object", GSL_EINVAL);
  77     }
  78
  79   for (i = 1; i < size; i++)
  80     {
  81       if (!(x_array[i-1] < x_array[i]))
  82         {
  83           GSL_ERROR ("x values must be monotonically increasing", GSL_EINVAL);
  84         }
  85     }
  86
  87   interp->xmin = x_array[0];
  88   interp->xmax = x_array[size - 1];
  89
  90   {
  91     int status = interp->type->init(interp->state, x_array, y_array, size);
  92     return status;
  93   }
  94 }
  95
  96 const char *
  97 gsl_interp_name(const gsl_interp * interp)
  98 {
  99   return interp->type->name;
 100 }
 101
 102 unsigned int
 103 gsl_interp_min_size(const gsl_interp * interp)
 104 {
 105   return interp->type->min_size;
 106 }
 107
 108 void
 109 gsl_interp_free (gsl_interp * interp)
 110 {
 111   if (interp->type->free)
 112     interp->type->free (interp->state);
 113   free (interp);
 114 }
 115
 116
 117
 118 int
 119 gsl_interp_eval_e (const gsl_interp * interp,
 120                    const double xa[], const double ya[], double x,
 121                    gsl_interp_accel * a, double *y)
 122 {
 123   if (x < interp->xmin)
 124     {
 125       *y = ya[0];
 126       return GSL_EDOM;
 127     }
 128   else if (x > interp->xmax)
 129     {
 130       *y = ya[interp->size - 1];
 131       return GSL_EDOM;
 132     }
 133
 134   return interp->type->eval (interp->state, xa, ya, interp->size, x, a, y);
 135 }
 136
 137 double
 138 gsl_interp_eval (const gsl_interp * interp,
 139                  const double xa[], const double ya[], double x,
 140                  gsl_interp_accel * a)
 141 {
 142   double y;
 143   int status = interp->type->eval (interp->state, xa, ya, interp->size, x, a, &y);
 144
 145   DISCARD_STATUS(status);
 146
 147   return y;
 148 }
 149
 150
 151 int
 152 gsl_interp_eval_deriv_e (const gsl_interp * interp,
 153                          const double xa[], const double ya[], double x,
 154                          gsl_interp_accel * a,
 155                          double *dydx)
 156 {
 157   if (x < interp->xmin)
 158     {
 159       *dydx = 0.0;
 160       return GSL_EDOM;
 161     }
 162   else if (x > interp->xmax)
 163     {
 164       *dydx = 0.0;
 165       return GSL_EDOM;
 166     }
 167
 168   return interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, dydx);
 169 }
 170
 171 double
 172 gsl_interp_eval_deriv (const gsl_interp * interp,
 173                        const double xa[], const double ya[], double x,
 174                        gsl_interp_accel * a)
 175 {
 176   double dydx;
 177   int status = interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, &dydx);
 178
 179   DISCARD_STATUS(status);
 180
 181   return dydx;
 182 }
 183
 184
 185 int
 186 gsl_interp_eval_deriv2_e (const gsl_interp * interp,
 187                           const double xa[], const double ya[], double x,
 188                           gsl_interp_accel * a,
 189                           double * d2)
 190 {
 191   if (x < interp->xmin)
 192     {
 193       *d2 = 0.0;
 194       return GSL_EDOM;
 195     }
 196   else if (x > interp->xmax)
 197     {
 198       *d2 = 0.0;
 199       return GSL_EDOM;
 200     }
 201
 202   return interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, d2);
 203 }
 204
 205 double
 206 gsl_interp_eval_deriv2 (const gsl_interp * interp,
 207                         const double xa[], const double ya[], double x,
 208                         gsl_interp_accel * a)
 209 {
 210   double d2;
 211   int status = interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, &d2);
 212
 213   DISCARD_STATUS(status);
 214
 215   return d2;
 216 }
 217
 218
 219 int
 220 gsl_interp_eval_integ_e (const gsl_interp * interp,
 221                          const double xa[], const double ya[],
 222                          double a, double b,
 223                          gsl_interp_accel * acc,
 224                          double * result)
 225 {
 226   if (a > b || a < interp->xmin || b > interp->xmax)
 227     {
 228       *result = 0.0;
 229       return GSL_EDOM;
 230     }
 231   else if(a == b)
 232     {
 233       *result = 0.0;
 234       return GSL_SUCCESS;
 235     }
 236
 237   return interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, result);
 238 }
 239
 240
 241 double
 242 gsl_interp_eval_integ (const gsl_interp * interp,
 243                        const double xa[], const double ya[],
 244                        double a, double b,
 245                        gsl_interp_accel * acc)
 246 {
 247   double result;
 248   int status = interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, &result);
 249
 250   DISCARD_STATUS(status);
 251
 252   return result;
 253 }
 254
 255