source: sasmodels/sasmodels/models/lib/jnd.c @ 3936ad3

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

Bessel functions from double-precison cephes has been added

  • Property mode set to 100644
File size: 1.8 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#define MACHEP 1.11022302462515654042E-16
51
52double jn( int n, double x )
53{
54double pkm2, pkm1, pk, xk, r, ans;
55int k, sign;
56
57if( n < 0 )
58        {
59        n = -n;
60        if( (n & 1) == 0 )      /* -1**n */
61                sign = 1;
62        else
63                sign = -1;
64        }
65else
66        sign = 1;
67
68if( x < 0.0 )
69        {
70        if( n & 1 )
71                sign = -sign;
72        x = -x;
73        }
74
75if( n == 0 )
76        return( sign * j0(x) );
77if( n == 1 )
78        return( sign * j1(x) );
79if( n == 2 )
80        return( sign * (2.0 * j1(x) / x  -  j0(x)) );
81
82if( x < MACHEP )
83        return( 0.0 );
84
85
86pk = 2 * (n + k);
87ans = pk;
88xk = x * x;
89
90do
91        {
92        pk -= 2.0;
93        ans = pk - (xk/ans);
94        }
95while( --k > 0 );
96ans = x/ans;
97
98/* backward recurrence */
99
100pk = 1.0;
101pkm1 = 1.0/ans;
102k = n-1;
103r = 2 * k;
104
105do
106        {
107        pkm2 = (pkm1 * r  -  pk * x) / x;
108        pk = pkm1;
109        pkm1 = pkm2;
110        r -= 2.0;
111        }
112while( --k > 0 );
113
114if( fabs(pk) > fabs(pkm1) )
115        ans = j1(x)/pk;
116else
117        ans = j0(x)/pkm1;
118return( sign * ans );
119}
120
Note: See TracBrowser for help on using the repository browser.