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