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 zero
7 
 
8 
	j0(x) returns the value of J0(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 J0 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 
	#5849 (19.22D)
23 
	#6549 (19.25D)
24 
	#6949 (19.41D)
25 
 
26 
	y0(x) returns the value of Y0(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, j0.
32 
 
33 
	The values of Y0 have not been checked
34 
	to more than ten places.
35 
 
36 
	Coefficients are from Hart & Cheney.
37 
	#6245 (18.78D)
38 
	#6549 (19.25D)
39 
	#6949 (19.41D)
40 
*/
41 
 
42 
static double pzero, qzero;
43 
static double tpi	= .6366197723675813430755350535e0;
44 
static double pio4	= .7853981633974483096156608458e0;
45 
static double p1[] = {
46 
	0.4933787251794133561816813446e21,
47 
	-.1179157629107610536038440800e21,
48 
	0.6382059341072356562289432465e19,
49 
	-.1367620353088171386865416609e18,
50 
	0.1434354939140344111664316553e16,
51 
	-.8085222034853793871199468171e13,
52 
	0.2507158285536881945555156435e11,
53 
	-.4050412371833132706360663322e8,
54 
	0.2685786856980014981415848441e5,
55 
};
56 
static double q1[] = {
57 
	0.4933787251794133562113278438e21,
58 
	0.5428918384092285160200195092e19,
59 
	0.3024635616709462698627330784e17,
60 
	0.1127756739679798507056031594e15,
61 
	0.3123043114941213172572469442e12,
62 
	0.6699987672982239671814028660e9,
63 
	0.1114636098462985378182402543e7,
64 
	0.1363063652328970604442810507e4,
65 
	1.0
66 
};
67 
static double p2[] = {
68 
	0.5393485083869438325262122897e7,
69 
	0.1233238476817638145232406055e8,
70 
	0.8413041456550439208464315611e7,
71 
	0.2016135283049983642487182349e7,
72 
	0.1539826532623911470917825993e6,
73 
	0.2485271928957404011288128951e4,
74 
	0.0,
75 
};
76 
static double q2[] = {
77 
	0.5393485083869438325560444960e7,
78 
	0.1233831022786324960844856182e8,
79 
	0.8426449050629797331554404810e7,
80 
	0.2025066801570134013891035236e7,
81 
	0.1560017276940030940592769933e6,
82 
	0.2615700736920839685159081813e4,
83 
	1.0,
84 
};
85 
static double p3[] = {
86 
	-.3984617357595222463506790588e4,
87 
	-.1038141698748464093880530341e5,
88 
	-.8239066313485606568803548860e4,
89 
	-.2365956170779108192723612816e4,
90 
	-.2262630641933704113967255053e3,
91 
	-.4887199395841261531199129300e1,
92 
	0.0,
93 
};
94 
static double q3[] = {
95 
	0.2550155108860942382983170882e6,
96 
	0.6667454239319826986004038103e6,
97 
	0.5332913634216897168722255057e6,
98 
	0.1560213206679291652539287109e6,
99 
	0.1570489191515395519392882766e5,
100 
	0.4087714673983499223402830260e3,
101 
	1.0,
102 
};
103 
static double p4[] = {
104 
	-.2750286678629109583701933175e20,
105 
	0.6587473275719554925999402049e20,
106 
	-.5247065581112764941297350814e19,
107 
	0.1375624316399344078571335453e18,
108 
	-.1648605817185729473122082537e16,
109 
	0.1025520859686394284509167421e14,
110 
	-.3436371222979040378171030138e11,
111 
	0.5915213465686889654273830069e8,
112 
	-.4137035497933148554125235152e5,
113 
};
114 
static double q4[] = {
115 
	0.3726458838986165881989980e21,
116 
	0.4192417043410839973904769661e19,
117 
	0.2392883043499781857439356652e17,
118 
	0.9162038034075185262489147968e14,
119 
	0.2613065755041081249568482092e12,
120 
	0.5795122640700729537480087915e9,
121 
	0.1001702641288906265666651753e7,
122 
	0.1282452772478993804176329391e4,
123 
	1.0,
124 
};
125 
 
126 
static void
127 
asympt(double arg)
128 
{
129 
	double zsq, n, d;
130 
	int i;
131 
 
132 
	zsq = 64 / (arg*arg);
133 
	for(n=0,d=0,i=6;i>=0;i--) {
134 
		n = n*zsq + p2[i];
135 
		d = d*zsq + q2[i];
136 
	}
137 
	pzero = n/d;
138 
	for(n=0,d=0,i=6;i>=0;i--) {
139 
		n = n*zsq + p3[i];
140 
		d = d*zsq + q3[i];
141 
	}
142 
	qzero = (8/arg)*(n/d);
143 
}
144 
 
145 
double
146 
j0(double arg)
147 
{
148 
	double argsq, n, d;
149 
	int i;
150 
 
151 
	if(arg < 0)
152 
		arg = -arg;
153 
	if(arg > 8) {
154 
		asympt(arg);
155 
		n = arg - pio4;
156 
		return sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n));
157 
	}
158 
	argsq = arg*arg;
159 
	for(n=0,d=0,i=8;i>=0;i--) {
160 
		n = n*argsq + p1[i];
161 
		d = d*argsq + q1[i];
162 
	}
163 
	return n/d;
164 
}
165 
 
166 
double
167 
y0(double arg)
168 
{
169 
	double argsq, n, d;
170 
	int i;
171 
 
172 
	errno = 0;
173 
	if(arg <= 0) {
174 
		errno = EDOM;
175 
		return(-HUGE_VAL);
176 
	}
177 
	if(arg > 8) {
178 
		asympt(arg);
179 
		n = arg - pio4;
180 
		return sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n));
181 
	}
182 
	argsq = arg*arg;
183 
	for(n=0,d=0,i=8;i>=0;i--) {
184 
		n = n*argsq + p4[i];
185 
		d = d*argsq + q4[i];
186 
	}
187 
	return n/d + tpi*j0(arg)*log(arg);
188 
}