1 #include <math.h>
2 #include <errno.h>
3
4 /*
5 floating point Bessel's function of
6 the first and second kinds and of
7 integer order.
8
9 int n;
10 double x;
11 jn(n,x);
12
13 returns the value of Jn(x) for all
14 integer values of n and all real values
15 of x.
16
17 There are no error returns.
18 Calls j0, j1.
19
20 For n=0, j0(x) is called,
21 for n=1, j1(x) is called,
22 for n<x, forward recursion us used starting
23 from values of j0(x) and j1(x).
24 for n>x, a continued fraction approximation to
25 j(n,x)/j(n-1,x) is evaluated and then backward
26 recursion is used starting from a supposed value
27 for j(n,x). The resulting value of j(0,x) is
28 compared with the actual value to correct the
29 supposed value of j(n,x).
30
31 yn(n,x) is similar in all respects, except
32 that forward recursion is used for all
33 values of n>1.
34 */
35
36 double j0(double);
37 double j1(double);
38 double y0(double);
39 double y1(double);
40
41 double
jn(int n,double x)42 jn(int n, double x)
43 {
44 int i;
45 double a, b, temp;
46 double xsq, t;
47
48 if(n < 0) {
49 n = -n;
50 x = -x;
51 }
52 if(n == 0)
53 return j0(x);
54 if(n == 1)
55 return j1(x);
56 if(x == 0)
57 return 0;
58 if(n > x)
59 goto recurs;
60
61 a = j0(x);
62 b = j1(x);
63 for(i=1; i<n; i++) {
64 temp = b;
65 b = (2*i/x)*b - a;
66 a = temp;
67 }
68 return b;
69
70 recurs:
71 xsq = x*x;
72 for(t=0,i=n+16; i>n; i--)
73 t = xsq/(2*i - t);
74 t = x/(2*n-t);
75
76 a = t;
77 b = 1;
78 for(i=n-1; i>0; i--) {
79 temp = b;
80 b = (2*i/x)*b - a;
81 a = temp;
82 }
83 return t*j0(x)/b;
84 }
85
86 double
yn(int n,double x)87 yn(int n, double x)
88 {
89 int i;
90 int sign;
91 double a, b, temp;
92
93 if (x <= 0) {
94 errno = EDOM;
95 return -HUGE_VAL;
96 }
97 sign = 1;
98 if(n < 0) {
99 n = -n;
100 if(n%2 == 1)
101 sign = -1;
102 }
103 if(n == 0)
104 return y0(x);
105 if(n == 1)
106 return sign*y1(x);
107
108 a = y0(x);
109 b = y1(x);
110 for(i=1; i<n; i++) {
111 temp = b;
112 b = (2*i/x)*b - a;
113 a = temp;
114 }
115 return sign*b;
116 }
117