git.gag.com Git - debian/gnuradio/blob - gnuradio-core/src/gen_interpolator_taps/simpson.c

   1 /* -*- c -*- */
   2 #include <math.h>
   3 #include <stdio.h>
   4
   5 #define EPS     (1.0e-5)
   6 #define JMAX    16
   7
   8 /*
   9  * Compute the Nth stage of refinement of an extended trapezoidal
  10  * rule.  FUNC is input as a pointer to a function to be integrated
  11  * between limits A and B.  When called with N = 1, the routine
  12  * returns the crudest estimate of the integral from A to B of f(x)
  13  * dx.  Subsequent calls with N=2,3,... (in that sequential order)
  14  * will improve the accuracy by adding 2**(N-2) additional interior
  15  * points.
  16  *
  17  * N.B., this function contains static state and IS NEITHER RENTRANT
  18  * NOR THREAD SAFE!
  19  */
  20
  21 double
  22 trapzd (double (*func)(double),
  23         double a, double b,
  24         int    n)
  25 {
  26   long double           x, tnm, sum, del;
  27   static long double    s;
  28   static int    it;
  29   int           j;
  30
  31   if (n == 1){
  32     it = 1;     /* # of points to add on the next call */
  33     s = 0.5 * (b - a) * (func(a) + func(b));
  34     return s;
  35   }
  36   else {
  37     tnm = it;
  38     del = (b-a)/tnm;            /* this is the spacing of the points to be added */
  39     x = a + 0.5*del;
  40     for (sum = 0.0, j = 1; j <= it; j++, x += del)
  41       sum += func(x);
  42     it *= 2;
  43     s = 0.5 * (s + (b-a) * sum/tnm);    /* replace s by it's refined value */
  44     return s;
  45   }
  46 }
  47
  48 /*
  49  * Returns the integral of the function FUNC from A to B.  The
  50  * parameters EPS can be set to the desired fractional accuracy and
  51  * JMAX so that 2**(JMAX-1) is the maximum allowed number of steps.
  52  * Integration is performed by Simpson's rule.
  53  */
  54
  55 double
  56 qsimp (double (*func)(double),
  57        double a,        /* lower limit */
  58        double b)        /* upper limit */
  59 {
  60   int           j;
  61   long double   s, st, ost, os;
  62
  63   ost = os = -1.0e30;
  64   for (j = 1; j <= JMAX; j++){
  65     st = trapzd (func, a, b, j);
  66     s = (4.0 * st - ost)/3.0;
  67     if (fabs (s - os) < EPS * fabs(os))
  68       return s;
  69     os = s;
  70     ost = st;
  71   }
  72   fprintf (stderr, "Too many steps in routine QSIMP\n");
  73   // exit (1);
  74   return s;
  75 }
  76