← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/special/cephes/jn.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/special/cephes/jn.c
 
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
 
45
Cephes Math Library Release 2.0:  April, 1987
 
46
Copyright 1984, 1987 by Stephen L. Moshier
 
47
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 
48
*/
 
49
#include "mconf.h"
 
50
#ifndef ANSIPROT
 
51
double fabs(), j0(), j1();
 
52
#endif
 
53
extern double MACHEP;
 
54
 
 
55
double jn( n, x )
 
56
int n;
 
57
double x;
 
58
{
 
59
double pkm2, pkm1, pk, xk, r, ans;
 
60
int k, sign;
 
61
 
 
62
if( n < 0 )
 
63
        {
 
64
        n = -n;
 
65
        if( (n & 1) == 0 )      /* -1**n */
 
66
                sign = 1;
 
67
        else
 
68
                sign = -1;
 
69
        }
 
70
else
 
71
        sign = 1;
 
72
 
 
73
if( x < 0.0 )
 
74
        {
 
75
        if( n & 1 )
 
76
                sign = -sign;
 
77
        x = -x;
 
78
        }
 
79
 
 
80
if( n == 0 )
 
81
        return( sign * j0(x) );
 
82
if( n == 1 )
 
83
        return( sign * j1(x) );
 
84
if( n == 2 )
 
85
        return( sign * (2.0 * j1(x) / x  -  j0(x)) );
 
86
 
 
87
if( x < MACHEP )
 
88
        return( 0.0 );
 
89
 
 
90
/* continued fraction */
 
91
#ifdef DEC
 
92
k = 56;
 
93
#else
 
94
k = 53;
 
95
#endif
 
96
 
 
97
pk = 2 * (n + k);
 
98
ans = pk;
 
99
xk = x * x;
 
100
 
 
101
do
 
102
        {
 
103
        pk -= 2.0;
 
104
        ans = pk - (xk/ans);
 
105
        }
 
106
while( --k > 0 );
 
107
ans = x/ans;
 
108
 
 
109
/* backward recurrence */
 
110
 
 
111
pk = 1.0;
 
112
pkm1 = 1.0/ans;
 
113
k = n-1;
 
114
r = 2 * k;
 
115
 
 
116
do
 
117
        {
 
118
        pkm2 = (pkm1 * r  -  pk * x) / x;
 
119
        pk = pkm1;
 
120
        pkm1 = pkm2;
 
121
        r -= 2.0;
 
122
        }
 
123
while( --k > 0 );
 
124
 
 
125
if( fabs(pk) > fabs(pkm1) )
 
126
        ans = j1(x)/pk;
 
127
else
 
128
        ans = j0(x)/pkm1;
 
129
return( sign * ans );
 
130
}