source: sasmodels/sasmodels/models/lib/jnd.c @ 0a9d219

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0a9d219 was 0a9d219, checked in by wojciech, 8 years ago

Single precisson cephes version has been added
  • Property mode set to 100644
File size: 2.7 KB
Line 
1/*                                                      jn.c
2 *
3 *      Bessel function of integer order
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * int n;
10 * double x, y, jn();
11 *
12 * y = jn( n, x );
13 *
14 *
15 *
16 * DESCRIPTION:
17 *
18 * Returns Bessel function of order n, where n is a
19 * (possibly negative) integer.
20 *
21 * The ratio of jn(x) to j0(x) is computed by backward
22 * recurrence.  First the ratio jn/jn-1 is found by a
23 * continued fraction expansion.  Then the recurrence
24 * relating successive orders is applied until j0 or j1 is
25 * reached.
26 *
27 * If n = 0 or 1 the routine for j0 or j1 is called
28 * directly.
29 *
30 *
31 *
32 * ACCURACY:
33 *
34 *                      Absolute error:
35 * arithmetic   range      # trials      peak         rms
36 *    DEC       0, 30        5500       6.9e-17     9.3e-18
37 *    IEEE      0, 30        5000       4.4e-16     7.9e-17
38 *
39 *
40 * Not suitable for large n or x. Use jv() instead.
41 *
42 */
43
44/*                                                      jn.c
45Cephes Math Library Release 2.8:  June, 2000
46Copyright 1984, 1987, 2000 by Stephen L. Moshier
47*/
48
49double jn( int n, double x );
50
51double jn( int n, double x ) {
52
53    const double MACHEP = 1.11022302462515654042E-16;
54    double pkm2, pkm1, pk, xk, r, ans, xinv;
55    int k, sign;
56
57    if( n < 0 ) {
58            n = -n;
59            if( (n & 1) == 0 )  /* -1**n */
60                    sign = 1;
61            else
62                    sign = -1;
63        }
64    else
65            sign = 1;
66
67    if( x < 0.0 ) {
68            if( n & 1 )
69                    sign = -sign;
70            x = -x;
71        }
72
73    if( n == 0 )
74            return( sign * j0(x) );
75    if( n == 1 )
76            return( sign * j1(x) );
77    if( n == 2 )
78            return( sign * (2.0 * j1(x) / x  -  j0(x)) );
79
80    if( x < MACHEP )
81            return( 0.0 );
82
83    if (FLOAT_SIZE > 4)
84        k = 53;
85    else
86        k = 24;
87
88    pk = 2 * (n + k);
89    ans = pk;
90    xk = x * x;
91
92    do {
93            pk -= 2.0;
94            ans = pk - (xk/ans);
95        } while( --k > 0 );
96
97    /* backward recurrence */
98
99    pk = 1.0;
100
101    if (FLOAT_SIZE > 4) {
102        ans = x/ans;
103        pkm1 = 1.0/ans;
104
105        k = n-1;
106        r = 2 * k;
107
108        do {
109            pkm2 = (pkm1 * r  -  pk * x) / x;
110                pk = pkm1;
111                pkm1 = pkm2;
112                r -= 2.0;
113            } while( --k > 0 );
114
115        if( fabs(pk) > fabs(pkm1) )
116                ans = j1(x)/pk;
117        else
118                ans = j0(x)/pkm1;
119
120            return( sign * ans );
121        }
122    else {
123        xinv = 1.0/x;
124        pkm1 = ans * xinv;
125        k = n-1;
126        r = (float )(2 * k);
127
128        do {
129                pkm2 = (pkm1 * r  -  pk * x) * xinv;
130                pk = pkm1;
131                pkm1 = pkm2;
132                r -= 2.0;
133            }
134        while( --k > 0 );
135
136        r = pk;
137        if( r < 0 )
138                r = -r;
139        ans = pkm1;
140        if( ans < 0 )
141                ans = -ans;
142
143        if( r > ans )  /* if( fabs(pk) > fabs(pkm1) ) */
144                ans = sign * j1(x)/pk;
145        else
146                ans = sign * j0(x)/pkm1;
147        return( ans );
148    }
149}
150
Note: See TracBrowser for help on using the repository browser.