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

   1 /* interpolation/akima.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 <math.h>
  25 #include <gsl/gsl_errno.h>
  26 #include "integ_eval.h"
  27 #include <gsl/gsl_interp.h>
  28
  29 typedef struct
  30 {
  31   double * b;
  32   double * c;
  33   double * d;
  34   double * _m;
  35 } akima_state_t;
  36
  37
  38 /* common creation */
  39 static void *
  40 akima_alloc (size_t size)
  41 {
  42   akima_state_t *state = (akima_state_t *) malloc (sizeof (akima_state_t));
  43
  44   if (state == NULL)
  45     {
  46       GSL_ERROR_NULL("failed to allocate space for state", GSL_ENOMEM);
  47     }
  48
  49   state->b = (double *) malloc (size * sizeof (double));
  50
  51   if (state->b == NULL)
  52     {
  53       free (state);
  54       GSL_ERROR_NULL("failed to allocate space for b", GSL_ENOMEM);
  55     }
  56
  57   state->c = (double *) malloc (size * sizeof (double));
  58
  59   if (state->c == NULL)
  60     {
  61       free (state->b);
  62       free (state);
  63       GSL_ERROR_NULL("failed to allocate space for c", GSL_ENOMEM);
  64     }
  65
  66   state->d = (double *) malloc (size * sizeof (double));
  67
  68   if (state->d == NULL)
  69     {
  70       free (state->c);
  71       free (state->b);
  72       free (state);
  73       GSL_ERROR_NULL("failed to allocate space for d", GSL_ENOMEM);
  74     }
  75
  76   state->_m = (double *) malloc ((size + 4) * sizeof (double));
  77
  78   if (state->_m == NULL)
  79     {
  80       free (state->d);
  81       free (state->c);
  82       free (state->b);
  83       free (state);
  84       GSL_ERROR_NULL("failed to allocate space for _m", GSL_ENOMEM);
  85     }
  86
  87   return state;
  88 }
  89
  90
  91 /* common calculation */
  92 static void
  93 akima_calc (const double x_array[], double b[],  double c[],  double d[], size_t size, double m[])
  94 {
  95   size_t i;
  96
  97   for (i = 0; i < (size - 1); i++)
  98     {
  99       const double NE = fabs (m[i + 1] - m[i]) + fabs (m[i - 1] - m[i - 2]);
 100       if (NE == 0.0)
 101         {
 102           b[i] = m[i];
 103           c[i] = 0.0;
 104           d[i] = 0.0;
 105         }
 106       else
 107         {
 108           const double h_i = x_array[i + 1] - x_array[i];
 109           const double NE_next = fabs (m[i + 2] - m[i + 1]) + fabs (m[i] - m[i - 1]);
 110           const double alpha_i = fabs (m[i - 1] - m[i - 2]) / NE;
 111           double alpha_ip1;
 112           double tL_ip1;
 113           if (NE_next == 0.0)
 114             {
 115               tL_ip1 = m[i];
 116             }
 117           else
 118             {
 119               alpha_ip1 = fabs (m[i] - m[i - 1]) / NE_next;
 120               tL_ip1 = (1.0 - alpha_ip1) * m[i] + alpha_ip1 * m[i + 1];
 121             }
 122           b[i] = (1.0 - alpha_i) * m[i - 1] + alpha_i * m[i];
 123           c[i] = (3.0 * m[i] - 2.0 * b[i] - tL_ip1) / h_i;
 124           d[i] = (b[i] + tL_ip1 - 2.0 * m[i]) / (h_i * h_i);
 125         }
 126     }
 127 }
 128
 129
 130 static int
 131 akima_init (void * vstate, const double x_array[], const double y_array[],
 132             size_t size)
 133 {
 134   akima_state_t *state = (akima_state_t *) vstate;
 135
 136   double * m = state->_m + 2; /* offset so we can address the -1,-2
 137                                  components */
 138
 139   size_t i;
 140   for (i = 0; i <= size - 2; i++)
 141     {
 142       m[i] = (y_array[i + 1] - y_array[i]) / (x_array[i + 1] - x_array[i]);
 143     }
 144
 145   /* non-periodic boundary conditions */
 146   m[-2] = 3.0 * m[0] - 2.0 * m[1];
 147   m[-1] = 2.0 * m[0] - m[1];
 148   m[size - 1] = 2.0 * m[size - 2] - m[size - 3];
 149   m[size] = 3.0 * m[size - 2] - 2.0 * m[size - 3];
 150
 151   akima_calc (x_array, state->b, state->c, state->d, size, m);
 152
 153   return GSL_SUCCESS;
 154 }
 155
 156
 157 static int
 158 akima_init_periodic (void * vstate,
 159                      const double x_array[],
 160                      const double y_array[],
 161                      size_t size)
 162 {
 163   akima_state_t *state = (akima_state_t *) vstate;
 164
 165   double * m = state->_m + 2; /* offset so we can address the -1,-2
 166                                  components */
 167
 168   size_t i;
 169   for (i = 0; i <= size - 2; i++)
 170     {
 171       m[i] = (y_array[i + 1] - y_array[i]) / (x_array[i + 1] - x_array[i]);
 172     }
 173
 174   /* periodic boundary conditions */
 175   m[-2] = m[size - 1 - 2];
 176   m[-1] = m[size - 1 - 1];
 177   m[size - 1] = m[0];
 178   m[size] = m[1];
 179
 180   akima_calc (x_array, state->b, state->c, state->d, size, m);
 181
 182   return GSL_SUCCESS;
 183 }
 184
 185 static void
 186 akima_free (void * vstate)
 187 {
 188   akima_state_t *state = (akima_state_t *) vstate;
 189
 190   free (state->b);
 191   free (state->c);
 192   free (state->d);
 193   free (state->_m);
 194   free (state);
 195 }
 196
 197
 198 static
 199 int
 200 akima_eval (const void * vstate,
 201             const double x_array[], const double y_array[], size_t size,
 202             double x,
 203             gsl_interp_accel * a,
 204             double *y)
 205 {
 206   const akima_state_t *state = (const akima_state_t *) vstate;
 207
 208   size_t index;
 209
 210   if (a != 0)
 211     {
 212       index = gsl_interp_accel_find (a, x_array, size, x);
 213     }
 214   else
 215     {
 216       index = gsl_interp_bsearch (x_array, x, 0, size - 1);
 217     }
 218
 219   /* evaluate */
 220   {
 221     const double x_lo = x_array[index];
 222     const double delx = x - x_lo;
 223     const double b = state->b[index];
 224     const double c = state->c[index];
 225     const double d = state->d[index];
 226     *y = y_array[index] + delx * (b + delx * (c + d * delx));
 227     return GSL_SUCCESS;
 228   }
 229 }
 230
 231
 232 static int
 233 akima_eval_deriv (const void * vstate,
 234                   const double x_array[], const double y_array[], size_t size,
 235                   double x,
 236                   gsl_interp_accel * a,
 237                   double *dydx)
 238 {
 239   const akima_state_t *state = (const akima_state_t *) vstate;
 240
 241   size_t index;
 242
 243   DISCARD_POINTER(y_array); /* prevent warning about unused parameter */
 244
 245   if (a != 0)
 246     {
 247       index = gsl_interp_accel_find (a, x_array, size, x);
 248     }
 249   else
 250     {
 251       index = gsl_interp_bsearch (x_array, x, 0, size - 1);
 252     }
 253
 254   /* evaluate */
 255   {
 256     double x_lo = x_array[index];
 257     double delx = x - x_lo;
 258     double b = state->b[index];
 259     double c = state->c[index];
 260     double d = state->d[index];
 261     *dydx = b + delx * (2.0 * c + 3.0 * d * delx);
 262     return GSL_SUCCESS;
 263   }
 264 }
 265
 266
 267 static
 268 int
 269 akima_eval_deriv2 (const void * vstate,
 270                    const double x_array[], const double y_array[], size_t size,
 271                    double x,
 272                    gsl_interp_accel * a,
 273                    double *y_pp)
 274 {
 275   const akima_state_t *state = (const akima_state_t *) vstate;
 276
 277   size_t index;
 278
 279   DISCARD_POINTER(y_array); /* prevent warning about unused parameter */
 280
 281   if (a != 0)
 282     {
 283       index = gsl_interp_accel_find (a, x_array, size, x);
 284     }
 285   else
 286     {
 287       index = gsl_interp_bsearch (x_array, x, 0, size - 1);
 288     }
 289
 290   /* evaluate */
 291   {
 292     const double x_lo = x_array[index];
 293     const double delx = x - x_lo;
 294     const double c = state->c[index];
 295     const double d = state->d[index];
 296     *y_pp = 2.0 * c + 6.0 * d * delx;
 297     return GSL_SUCCESS;
 298   }
 299 }
 300
 301
 302 static
 303 int
 304 akima_eval_integ (const void * vstate,
 305                   const double x_array[], const double y_array[], size_t size,
 306                   gsl_interp_accel * acc,
 307                   double a, double b,
 308                   double * result)
 309 {
 310   const akima_state_t *state = (const akima_state_t *) vstate;
 311
 312   size_t i, index_a, index_b;
 313
 314   if (acc != 0)
 315     {
 316       index_a = gsl_interp_accel_find (acc, x_array, size, a);
 317       index_b = gsl_interp_accel_find (acc, x_array, size, b);
 318     }
 319   else
 320     {
 321       index_a = gsl_interp_bsearch (x_array, a, 0, size - 1);
 322       index_b = gsl_interp_bsearch (x_array, b, 0, size - 1);
 323     }
 324
 325   *result = 0.0;
 326
 327   /* interior intervals */
 328
 329   for(i=index_a; i<=index_b; i++) {
 330     const double x_hi = x_array[i + 1];
 331     const double x_lo = x_array[i];
 332     const double y_lo = y_array[i];
 333     const double dx = x_hi - x_lo;
 334     if(dx != 0.0) {
 335
 336       if (i == index_a || i == index_b)
 337         {
 338           double x1 = (i == index_a) ? a : x_lo;
 339           double x2 = (i == index_b) ? b : x_hi;
 340           *result += integ_eval (y_lo, state->b[i], state->c[i], state->d[i],
 341                                  x_lo, x1, x2);
 342         }
 343       else
 344         {
 345           *result += dx * (y_lo
 346                            + dx*(0.5*state->b[i]
 347                                  + dx*(state->c[i]/3.0
 348                                        + 0.25*state->d[i]*dx)));
 349         }
 350     }
 351     else {
 352       *result = 0.0;
 353       return GSL_FAILURE;
 354     }
 355   }
 356
 357   return GSL_SUCCESS;
 358 }
 359
 360
 361 static const gsl_interp_type akima_type =
 362 {
 363   "akima",
 364   5,
 365   &akima_alloc,
 366   &akima_init,
 367   &akima_eval,
 368   &akima_eval_deriv,
 369   &akima_eval_deriv2,
 370   &akima_eval_integ,
 371   &akima_free
 372 };
 373
 374 const gsl_interp_type * gsl_interp_akima = &akima_type;
 375
 376 static const gsl_interp_type akima_periodic_type =
 377 {
 378   "akima-periodic",
 379   5,
 380   &akima_alloc,
 381   &akima_init_periodic,
 382   &akima_eval,
 383   &akima_eval_deriv,
 384   &akima_eval_deriv2,
 385   &akima_eval_integ,
 386   &akima_free
 387 };
 388
 389 const gsl_interp_type * gsl_interp_akima_periodic = &akima_periodic_type;