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

#### Defined functions

jn defined in line 44; never used
yn defined in line 86; never used

#### Defined variables

sccsid defined in line 2; never used
zero defined in line 41; used 2 times
• in line 98(2)
 Last modified: 1985-08-21 Generated: 2016-12-26 Generated by src2html V0.67 page hit count: 992