1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
|
/* specfunc/bessel_sequence.c
*
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/* Author: G. Jungman */
#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_bessel.h>
#define DYDX_p(p,u,x) (-(p)/(x) + (((nu)*(nu))/((x)*(x))-1.0)*(u))
#define DYDX_u(p,u,x) (p)
static int
rk_step(double nu, double x, double dx, double * Jp, double * J)
{
double p_0 = *Jp;
double u_0 = *J;
double p_1 = dx * DYDX_p(p_0, u_0, x);
double u_1 = dx * DYDX_u(p_0, u_0, x);
double p_2 = dx * DYDX_p(p_0 + 0.5*p_1, u_0 + 0.5*u_1, x + 0.5*dx);
double u_2 = dx * DYDX_u(p_0 + 0.5*p_1, u_0 + 0.5*u_1, x + 0.5*dx);
double p_3 = dx * DYDX_p(p_0 + 0.5*p_2, u_0 + 0.5*u_2, x + 0.5*dx);
double u_3 = dx * DYDX_u(p_0 + 0.5*p_2, u_0 + 0.5*u_2, x + 0.5*dx);
double p_4 = dx * DYDX_p(p_0 + p_3, u_0 + u_3, x + dx);
double u_4 = dx * DYDX_u(p_0 + p_3, u_0 + u_3, x + dx);
*Jp = p_0 + p_1/6.0 + p_2/3.0 + p_3/3.0 + p_4/6.0;
*J = u_0 + u_1/6.0 + u_2/3.0 + u_3/3.0 + u_4/6.0;
return GSL_SUCCESS;
}
int
gsl_sf_bessel_sequence_Jnu_e(double nu, gsl_mode_t mode, size_t size, double * v)
{
/* CHECK_POINTER(v) */
if(nu < 0.0) {
GSL_ERROR ("domain error", GSL_EDOM);
}
else if(size == 0) {
GSL_ERROR ("error", GSL_EINVAL);
}
else {
const gsl_prec_t goal = GSL_MODE_PREC(mode);
const double dx_array[] = { 0.001, 0.03, 0.1 }; /* double, single, approx */
const double dx_nominal = dx_array[goal];
const int cnu = (int) ceil(nu);
const double nu13 = pow(nu,1.0/3.0);
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 };
const double x_small = ( nu >= 10.0 ? nu - nu13 : smalls[cnu] );
gsl_sf_result J0, J1;
double Jp, J;
double x;
size_t i = 0;
/* Calculate the first point. */
x = v[0];
gsl_sf_bessel_Jnu_e(nu, x, &J0);
v[0] = J0.val;
++i;
/* Step over the idiot case where the
* first point was actually zero.
*/
if(x == 0.0) {
if(v[1] <= x) {
/* Strict ordering failure. */
GSL_ERROR ("error", GSL_EFAILED);
}
x = v[1];
gsl_sf_bessel_Jnu_e(nu, x, &J0);
v[1] = J0.val;
++i;
}
/* Calculate directly as long as the argument
* is small. This is necessary because the
* integration is not very good there.
*/
while(v[i] < x_small && i < size) {
if(v[i] <= x) {
/* Strict ordering failure. */
GSL_ERROR ("error", GSL_EFAILED);
}
x = v[i];
gsl_sf_bessel_Jnu_e(nu, x, &J0);
v[i] = J0.val;
++i;
}
/* At this point we are ready to integrate.
* The value of x is the last calculated
* point, which has the value J0; v[i] is
* the next point we need to calculate. We
* calculate nu+1 at x as well to get
* the derivative, then we go forward.
*/
gsl_sf_bessel_Jnu_e(nu+1.0, x, &J1);
J = J0.val;
Jp = -J1.val + nu/x * J0.val;
while(i < size) {
const double dv = v[i] - x;
const int Nd = (int) ceil(dv/dx_nominal);
const double dx = dv / Nd;
double xj;
int j;
if(v[i] <= x) {
/* Strict ordering failure. */
GSL_ERROR ("error", GSL_EFAILED);
}
/* Integrate over interval up to next sample point.
*/
for(j=0, xj=x; j<Nd; j++, xj += dx) {
rk_step(nu, xj, dx, &Jp, &J);
}
/* Go to next interval. */
x = v[i];
v[i] = J;
++i;
}
return GSL_SUCCESS;
}
}
|
