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

#### Defined variables

errno defined in line 38; used 1 times
• in line 90
 Last modified: 1985-06-05 Generated: 2016-12-26 Generated by src2html V0.67 page hit count: 1880  