← Back to branch summary

~ubuntu-branches/ubuntu/lucid/pdl/lucid

« back to all changes in this revision

Viewing changes to Basic/Math/jn.c

  • Committer: Bazaar Package Importer
  • Author(s): Ben Gertzfield
  • Date: 2002-04-08 18:47:16 UTC
  • Revision ID: james.westby@ubuntu.com-20020408184716-0hf64dc96kin3htp
Tags: upstream-2.3.2
Import upstream version 2.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Basic/Math/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
 
 
51
extern double MACHEP;
 
52
 
 
53
double jn( n, x )
 
54
int n;
 
55
double x;
 
56
{
 
57
double pkm2, pkm1, pk, xk, r, ans;
 
58
int k, sign;
 
59
double fabs(), j0(), j1();
 
60
 
 
61
if( n < 0 )
 
62
        {
 
63
        n = -n;
 
64
        if( (n & 1) == 0 )      /* -1**n */
 
65
                sign = 1;
 
66
        else
 
67
                sign = -1;
 
68
        }
 
69
else
 
70
        sign = 1;
 
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
/* continued fraction */
 
84
#ifdef DEC
 
85
k = 56;
 
86
#else
 
87
k = 53;
 
88
#endif
 
89
 
 
90
pk = 2 * (n + k);
 
91
ans = pk;
 
92
xk = x * x;
 
93
 
 
94
do
 
95
        {
 
96
        pk -= 2.0;
 
97
        ans = pk - (xk/ans);
 
98
        }
 
99
while( --k > 0 );
 
100
ans = x/ans;
 
101
 
 
102
/* backward recurrence */
 
103
 
 
104
pk = 1.0;
 
105
pkm1 = 1.0/ans;
 
106
k = n-1;
 
107
r = 2 * k;
 
108
 
 
109
do
 
110
        {
 
111
        pkm2 = (pkm1 * r  -  pk * x) / x;
 
112
        pk = pkm1;
 
113
        pkm1 = pkm2;
 
114
        r -= 2.0;
 
115
        }
 
116
while( --k > 0 );
 
117
 
 
118
if( fabs(pk) > fabs(pkm1) )
 
119
        ans = j1(x)/pk;
 
120
else
 
121
        ans = j0(x)/pkm1;
 
122
return( sign * ans );
 
123
}