Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 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
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
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 
}