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