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