source: sasmodels/sasmodels/models/lib/sas_JN.c @ 4937980

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 4937980 was 4937980, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

Merge remote-tracking branch 'origin/master' into polydisp

Conflicts:

sasmodels/models/flexible_cylinder.c
  • 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 sas_JN( int n, double x );
50
51double sas_JN( int n, double x ) {
52
53    const double MACHEP = 1.11022302462515654042E-16;
54    double pkm2, pkm1, pk, xk, r, ans;
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 * sas_J0(x) );
75    if( n == 1 )
76            return( sign * sas_J1(x) );
77    if( n == 2 )
78            return( sign * (2.0 * sas_J1(x) / x  -  sas_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    #endif
88
89    pk = 2 * (n + k);
90    ans = pk;
91    xk = x * x;
92
93    do {
94            pk -= 2.0;
95            ans = pk - (xk/ans);
96        } while( --k > 0 );
97
98    /* backward recurrence */
99
100    pk = 1.0;
101
102    #if FLOAT_SIZE > 4
103        ans = x/ans;
104        pkm1 = 1.0/ans;
105
106        k = n-1;
107        r = 2 * k;
108
109        do {
110            pkm2 = (pkm1 * r  -  pk * x) / x;
111                pk = pkm1;
112                pkm1 = pkm2;
113                r -= 2.0;
114            } while( --k > 0 );
115
116        if( fabs(pk) > fabs(pkm1) )
117                ans = sas_J1(x)/pk;
118        else
119                ans = sas_J0(x)/pkm1;
120
121            return( sign * ans );
122
123    #else
124        const double xinv = 1.0/x;
125        pkm1 = ans * xinv;
126        k = n-1;
127        r = (float )(2 * k);
128
129        do {
130                pkm2 = (pkm1 * r  -  pk * x) * xinv;
131                pk = pkm1;
132                pkm1 = pkm2;
133                r -= 2.0;
134            }
135        while( --k > 0 );
136
137        r = pk;
138        if( r < 0 )
139                r = -r;
140        ans = pkm1;
141        if( ans < 0 )
142                ans = -ans;
143
144        if( r > ans )  /* if( fabs(pk) > fabs(pkm1) ) */
145                ans = sign * sas_J1(x)/pk;
146        else
147                ans = sign * sas_J0(x)/pkm1;
148        return( ans );
149    #endif
150}
151
Note: See TracBrowser for help on using the repository browser.