htswanalysis/MACS/lib/gsl/gsl-1.11/cheb/init.c

   1 /* cheb/init.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 #include <config.h>
  21 #include <stdlib.h>
  22 #include <gsl/gsl_math.h>
  23 #include <gsl/gsl_errno.h>
  24 #include <gsl/gsl_chebyshev.h>
  25
  26 /*-*-*-*-*-*-*-*-*-*-*-* Allocators *-*-*-*-*-*-*-*-*-*-*-*/
  27
  28 gsl_cheb_series *
  29 gsl_cheb_alloc(const size_t order)
  30 {
  31   gsl_cheb_series * cs = (gsl_cheb_series *) malloc(sizeof(gsl_cheb_series));
  32
  33   if(cs == 0) {
  34     GSL_ERROR_VAL("failed to allocate gsl_cheb_series struct", GSL_ENOMEM, 0);
  35   }
  36
  37   cs->order    = order;
  38   cs->order_sp = order;
  39
  40   cs->c = (double *) malloc((order+1) * sizeof(double));
  41
  42   if(cs->c == 0) {
  43     GSL_ERROR_VAL("failed to allocate cheb coefficients", GSL_ENOMEM, 0);
  44   }
  45
  46   cs->f = (double *) malloc((order+1) * sizeof(double));
  47
  48   if(cs->f == 0) {
  49     GSL_ERROR_VAL("failed to allocate cheb function space", GSL_ENOMEM, 0);
  50   }
  51
  52   return cs;
  53 }
  54
  55
  56 void gsl_cheb_free(gsl_cheb_series * cs)
  57 {
  58   free(cs->f);
  59   free(cs->c);
  60   free(cs);
  61 }
  62
  63 /*-*-*-*-*-*-*-*-*-*-*-* Initializer *-*-*-*-*-*-*-*-*-*-*-*/
  64
  65 int gsl_cheb_init(gsl_cheb_series * cs, const gsl_function *func,
  66                   const double a, const double b)
  67 {
  68   size_t k, j;
  69
  70   if(a >= b) {
  71     GSL_ERROR_VAL("null function interval [a,b]", GSL_EDOM, 0);
  72   }
  73   cs->a = a;
  74   cs->b = b;
  75   /* cs->err = 0.0; */
  76
  77   {
  78     double bma = 0.5 * (cs->b - cs->a);
  79     double bpa = 0.5 * (cs->b + cs->a);
  80     double fac = 2.0/(cs->order +1.0);
  81
  82     for(k = 0; k<=cs->order; k++) {
  83       double y = cos(M_PI * (k+0.5)/(cs->order+1));
  84       cs->f[k] = GSL_FN_EVAL(func, (y*bma + bpa));
  85     }
  86
  87     for(j = 0; j<=cs->order; j++) {
  88       double sum = 0.0;
  89       for(k = 0; k<=cs->order; k++)
  90         sum += cs->f[k]*cos(M_PI * j*(k+0.5)/(cs->order+1));
  91       cs->c[j] = fac * sum;
  92     }
  93
  94   }
  95   return GSL_SUCCESS;
  96 }
  97
  98
  99
 100
 101