htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/bessel_sequence.c

   1 /* specfunc/bessel_sequence.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 <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_bessel.h>
  26
  27
  28 #define DYDX_p(p,u,x) (-(p)/(x) + (((nu)*(nu))/((x)*(x))-1.0)*(u))
  29 #define DYDX_u(p,u,x) (p)
  30
  31 static int
  32 rk_step(double nu, double x, double dx, double * Jp, double * J)
  33 {
  34   double p_0 = *Jp;
  35   double u_0 = *J;
  36
  37   double p_1 = dx * DYDX_p(p_0, u_0, x);
  38   double u_1 = dx * DYDX_u(p_0, u_0, x);
  39
  40   double p_2 = dx * DYDX_p(p_0 + 0.5*p_1, u_0 + 0.5*u_1, x + 0.5*dx);
  41   double u_2 = dx * DYDX_u(p_0 + 0.5*p_1, u_0 + 0.5*u_1, x + 0.5*dx);
  42
  43   double p_3 = dx * DYDX_p(p_0 + 0.5*p_2, u_0 + 0.5*u_2, x + 0.5*dx);
  44   double u_3 = dx * DYDX_u(p_0 + 0.5*p_2, u_0 + 0.5*u_2, x + 0.5*dx);
  45
  46   double p_4 = dx * DYDX_p(p_0 + p_3, u_0 + u_3, x + dx);
  47   double u_4 = dx * DYDX_u(p_0 + p_3, u_0 + u_3, x + dx);
  48
  49   *Jp = p_0 + p_1/6.0 + p_2/3.0 + p_3/3.0 + p_4/6.0;
  50   *J  = u_0 + u_1/6.0 + u_2/3.0 + u_3/3.0 + u_4/6.0;
  51
  52   return GSL_SUCCESS;
  53 }
  54
  55
  56 int
  57 gsl_sf_bessel_sequence_Jnu_e(double nu, gsl_mode_t mode, size_t size, double * v)
  58 {
  59   /* CHECK_POINTER(v) */
  60
  61   if(nu < 0.0) {
  62     GSL_ERROR ("domain error", GSL_EDOM);
  63   }
  64   else if(size == 0) {
  65     GSL_ERROR ("error", GSL_EINVAL);
  66   }
  67   else {
  68     const gsl_prec_t goal   = GSL_MODE_PREC(mode);
  69     const double dx_array[] = { 0.001, 0.03, 0.1 }; /* double, single, approx */
  70     const double dx_nominal = dx_array[goal];
  71
  72     const int cnu = (int) ceil(nu);
  73     const double nu13 = pow(nu,1.0/3.0);
  74     const double smalls[] = { 0.01, 0.02, 0.4, 0.7, 1.3, 2.0, 2.5, 3.2, 3.5, 4.5, 6.0 };
  75     const double x_small = ( nu >= 10.0 ? nu - nu13 : smalls[cnu] );
  76
  77     gsl_sf_result J0, J1;
  78     double Jp, J;
  79     double x;
  80     size_t i = 0;
  81
  82     /* Calculate the first point. */
  83     x = v[0];
  84     gsl_sf_bessel_Jnu_e(nu, x, &J0);
  85     v[0] = J0.val;
  86     ++i;
  87
  88     /* Step over the idiot case where the
  89      * first point was actually zero.
  90      */
  91     if(x == 0.0) {
  92       if(v[1] <= x) {
  93         /* Strict ordering failure. */
  94         GSL_ERROR ("error", GSL_EFAILED);
  95       }
  96       x = v[1];
  97       gsl_sf_bessel_Jnu_e(nu, x, &J0);
  98       v[1] = J0.val;
  99       ++i;
 100     }
 101
 102     /* Calculate directly as long as the argument
 103      * is small. This is necessary because the
 104      * integration is not very good there.
 105      */
 106     while(v[i] < x_small && i < size) {
 107       if(v[i] <= x) {
 108         /* Strict ordering failure. */
 109         GSL_ERROR ("error", GSL_EFAILED);
 110       }
 111       x = v[i];
 112       gsl_sf_bessel_Jnu_e(nu, x, &J0);
 113       v[i] = J0.val;
 114       ++i;
 115     }
 116
 117     /* At this point we are ready to integrate.
 118      * The value of x is the last calculated
 119      * point, which has the value J0; v[i] is
 120      * the next point we need to calculate. We
 121      * calculate nu+1 at x as well to get
 122      * the derivative, then we go forward.
 123      */
 124     gsl_sf_bessel_Jnu_e(nu+1.0, x, &J1);
 125     J  = J0.val;
 126     Jp = -J1.val + nu/x * J0.val;
 127
 128     while(i < size) {
 129       const double dv = v[i] - x;
 130       const int Nd    = (int) ceil(dv/dx_nominal);
 131       const double dx = dv / Nd;
 132       double xj;
 133       int j;
 134
 135       if(v[i] <= x) {
 136         /* Strict ordering failure. */
 137         GSL_ERROR ("error", GSL_EFAILED);
 138       }
 139
 140       /* Integrate over interval up to next sample point.
 141        */
 142       for(j=0, xj=x; j<Nd; j++, xj += dx) {
 143         rk_step(nu, xj, dx, &Jp, &J);
 144       }
 145
 146       /* Go to next interval. */
 147       x = v[i];
 148       v[i] = J;
 149       ++i;
 150     }
 151
 152     return GSL_SUCCESS;
 153   }
 154 }