home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / interpolation / linear.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  4.6 KB  |  199 lines
  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 2 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., 675 Mass Ave, Cambridge, MA 02139, 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;
  199.