Date: Sun, 12 Feb 95 18:49:29 +0000

	This is a transcription of algorithm #5 of Collected Algorithms from
ACM. I took a few liberties to express it in modern ASCII codes, namely,
the 'not equal' construct is expressed as =/= and the exponentiation
is expressed as a ^ character. Hence ((X/2) ^ (n+2*s)) means
             n+2*s
	X / 2

multiplication operators have been substituted by an * as in most modern
languages.

			Jose R. Valverde
		European Bioinformatics Institute

			txomsy@ebi.ac.uk

------------------------------------------------------------------------

5. BESSEL FUNCTION I, SERIES EXPANSION
   Dorothea S. Clarke
   General Electric Co., FPLD, Cincinnati 15, Ohio

comment	Compute the Bessel function In(X) when n and X
	are within the bounds of the series expansion.
	The procedure calling statement gives n, X and an
	absolute tolerance delta for determining the point at
	which the terms of the summation become insig-
	nificant. Special case: I0(0) = 1;

procedure	I(n, X, delta) =:(Is)
begin
I:		s := 0 ; sum := 0
if		(n =/= 0) ; go to STRT
if		(X = 0) ; begin Is := 1 ; return end
		summ := 1 ; go to SURE
STRT:		sfac := 1
if		(s = 0) ; go to HRE
for		t := 1 (1) s
		sfac := sfac * t
HRE:		snfac := sfac
for		t := s + 1 (1) s + n
		snfac := snfac * t
		summ := sum + ((X/2) ^ (n+2*s)) / sfac * snfac)
SURE: if	(delta < abs(summ - sum))
begin		s := s + 1 ; sum := summ ; go to STRT end
		Is := summ ; return
end