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 
	floating point Bessel's function
5 
	of the first and second kinds
6 
	of order one
7 
 
8 
	j1(x) returns the value of J1(x)
9 
	for all real values of x.
10 
 
11 
	There are no error returns.
12 
	Calls sin, cos, sqrt.
13 
 
14 
	There is a niggling bug in J1 which
15 
	causes errors up to 2e-16 for x in the
16 
	interval [-8,8].
17 
	The bug is caused by an inappropriate order
18 
	of summation of the series.  rhm will fix it
19 
	someday.
20 
 
21 
	Coefficients are from Hart & Cheney.
22 
	#6050 (20.98D)
23 
	#6750 (19.19D)
24 
	#7150 (19.35D)
25 
 
26 
	y1(x) returns the value of Y1(x)
27 
	for positive real values of x.
28 
	For x<=0, error number EDOM is set and a
29 
	large negative value is returned.
30 
 
31 
	Calls sin, cos, sqrt, log, j1.
32 
 
33 
	The values of Y1 have not been checked
34 
	to more than ten places.
35 
 
36 
	Coefficients are from Hart & Cheney.
37 
	#6447 (22.18D)
38 
	#6750 (19.19D)
39 
	#7150 (19.35D)
40 
*/
41 
 
42 
static double pzero, qzero;
43 
static double tpi	= .6366197723675813430755350535e0;
44 
static double pio4	= .7853981633974483096156608458e0;
45 
static double p1[] = {
46 
	0.581199354001606143928050809e21,
47 
	-.6672106568924916298020941484e20,
48 
	0.2316433580634002297931815435e19,
49 
	-.3588817569910106050743641413e17,
50 
	0.2908795263834775409737601689e15,
51 
	-.1322983480332126453125473247e13,
52 
	0.3413234182301700539091292655e10,
53 
	-.4695753530642995859767162166e7,
54 
	0.2701122710892323414856790990e4,
55 
};
56 
static double q1[] = {
57 
	0.1162398708003212287858529400e22,
58 
	0.1185770712190320999837113348e20,
59 
	0.6092061398917521746105196863e17,
60 
	0.2081661221307607351240184229e15,
61 
	0.5243710262167649715406728642e12,
62 
	0.1013863514358673989967045588e10,
63 
	0.1501793594998585505921097578e7,
64 
	0.1606931573481487801970916749e4,
65 
	1.0,
66 
};
67 
static double p2[] = {
68 
	-.4435757816794127857114720794e7,
69 
	-.9942246505077641195658377899e7,
70 
	-.6603373248364939109255245434e7,
71 
	-.1523529351181137383255105722e7,
72 
	-.1098240554345934672737413139e6,
73 
	-.1611616644324610116477412898e4,
74 
	0.0,
75 
};
76 
static double q2[] = {
77 
	-.4435757816794127856828016962e7,
78 
	-.9934124389934585658967556309e7,
79 
	-.6585339479723087072826915069e7,
80 
	-.1511809506634160881644546358e7,
81 
	-.1072638599110382011903063867e6,
82 
	-.1455009440190496182453565068e4,
83 
	1.0,
84 
};
85 
static double p3[] = {
86 
	0.3322091340985722351859704442e5,
87 
	0.8514516067533570196555001171e5,
88 
	0.6617883658127083517939992166e5,
89 
	0.1849426287322386679652009819e5,
90 
	0.1706375429020768002061283546e4,
91 
	0.3526513384663603218592175580e2,
92 
	0.0,
93 
};
94 
static double q3[] = {
95 
	0.7087128194102874357377502472e6,
96 
	0.1819458042243997298924553839e7,
97 
	0.1419460669603720892855755253e7,
98 
	0.4002944358226697511708610813e6,
99 
	0.3789022974577220264142952256e5,
100 
	0.8638367769604990967475517183e3,
101 
	1.0,
102 
};
103 
static double p4[] = {
104 
	-.9963753424306922225996744354e23,
105 
	0.2655473831434854326894248968e23,
106 
	-.1212297555414509577913561535e22,
107 
	0.2193107339917797592111427556e20,
108 
	-.1965887462722140658820322248e18,
109 
	0.9569930239921683481121552788e15,
110 
	-.2580681702194450950541426399e13,
111 
	0.3639488548124002058278999428e10,
112 
	-.2108847540133123652824139923e7,
113 
	0.0,
114 
};
115 
static double q4[] = {
116 
	0.5082067366941243245314424152e24,
117 
	0.5435310377188854170800653097e22,
118 
	0.2954987935897148674290758119e20,
119 
	0.1082258259408819552553850180e18,
120 
	0.2976632125647276729292742282e15,
121 
	0.6465340881265275571961681500e12,
122 
	0.1128686837169442121732366891e10,
123 
	0.1563282754899580604737366452e7,
124 
	0.1612361029677000859332072312e4,
125 
	1.0,
126 
};
127 
 
128 
static void
129 
asympt(double arg)
130 
{
131 
	double zsq, n, d;
132 
	int i;
133 
 
134 
	zsq = 64/(arg*arg);
135 
	for(n=0,d=0,i=6;i>=0;i--) {
136 
		n = n*zsq + p2[i];
137 
		d = d*zsq + q2[i];
138 
	}
139 
	pzero = n/d;
140 
	for(n=0,d=0,i=6;i>=0;i--) {
141 
		n = n*zsq + p3[i];
142 
		d = d*zsq + q3[i];
143 
	}
144 
	qzero = (8/arg)*(n/d);
145 
}
146 
 
147 
double
148 
j1(double arg)
149 
{
150 
	double xsq, n, d, x;
151 
	int i;
152 
 
153 
	x = arg;
154 
	if(x < 0)
155 
		x = -x;
156 
	if(x > 8) {
157 
		asympt(x);
158 
		n = x - 3*pio4;
159 
		n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
160 
		if(arg < 0)
161 
			n = -n;
162 
		return n;
163 
	}
164 
	xsq = x*x;
165 
	for(n=0,d=0,i=8;i>=0;i--) {
166 
		n = n*xsq + p1[i];
167 
		d = d*xsq + q1[i];
168 
	}
169 
	return arg*n/d;
170 
}
171 
 
172 
double
173 
y1(double arg)
174 
{
175 
	double xsq, n, d, x;
176 
	int i;
177 
 
178 
	errno = 0;
179 
	x = arg;
180 
	if(x <= 0) {
181 
		errno = EDOM;
182 
		return -HUGE_VAL;
183 
	}
184 
	if(x > 8) {
185 
		asympt(x);
186 
		n = x - 3*pio4;
187 
		return sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n));
188 
	}
189 
	xsq = x*x;
190 
	for(n=0,d=0,i=9;i>=0;i--) {
191 
		n = n*xsq + p4[i];
192 
		d = d*xsq + q4[i];
193 
	}
194 
	return x*n/d + tpi*(j1(x)*log(x)-1/x);
195 
}