/************************************************************************/
/* */
/* ALGORITHM 21 - BESSEL FUNCTION FOR A SET OF INTEGER ORDERS */
/* */
/* W. Borsch-Supan, National Bureau of Standards, Washington 25, D.C. */
/* */
/* Modified 16 November 1964, by J. Stafford, Westland Aircraft Ltd, */
/* Saunders-Roe Division, East Cowes, Isle of Wight, England. */
/* */
/************************************************************************/
/* Translated to PL/I by R. A. Vowels, 2 July 2006. */
/* Works 5 July 2006. */
/* This procedure computes the values of the Bessel */
/* functions Jp(x) for real argument x and the set of all integer */
/* orders from 0 up to n and stores these values into the array J, */
/* whose subscript bounds should include the integers from 0 up */
/* to n. n must be non-negative. */
/* The computation is done by applying the recurrence formula */
/* backward from p = k down to p = 0 as described in Mathematical */
/* Tables Aids Comp. 11 (1957), 255-257. k is chosen to yield */
/* errors less than 1e-5 */
/* approximately after the first application of the recurrence. The */
/* recurrence is repeated with a larger k until the difference between */
/* the results of the two last recurrences doesn't exceed the */
/* given bound eps > 0. The steps in increasing k are chosen in */
/* such a way that the errors decrease at least by a factor of */
/* approximately 1e-5. Protection against overflow - essential with */
/* the method - is provided in an alteration by J. Stafford. */
BESSELSETINT:
PROCEDURE (x, n, eps, J);
DECLARE (x, eps, J(*)) FLOAT (18);
DECLARE n FIXED BINARY;
DECLARE (dist, rec0, rec1, rec2, sum, max, err, Threshold,
Jbar(HBOUND(J,1)) ) FLOAT (18);
DECLARE (k, m, p) FIXED BINARY;
DECLARE s BIT(1);
IF LBOUND(J,1) ^= 0 THEN
DO;
PUT SKIP LIST ('***ERROR, lower bound of J must be zero'); STOP;
END;
IF x = 0 THEN
DO;
J(0) = 1;
DO p = 1 TO n;
J(p) = 0;
END;
RETURN;
END;
IF ABS(x) >= 8.0E0 THEN
dist = 5*ABS(x)**(1.00000E0/3);
ELSE
dist = 10;
IF ABS(x) >= n THEN
k = ABS(x) + dist + 1;
ELSE
k = n + dist + 1;
s = '0'b;
/* Stage 1: Calculate trial J values. */
/* For the recurrence formula, two start values are required, */
/* namely 0 and 1 (stored in rec0 and rec1 respectively). */
/* the recurrence formula is J(p-1) = 2*p/x * J(p) - J(p+1) */
/* [here, J(p) denotes Jp(x)]. */
/* In the following loop, J0 + 2J2 + 2J4 + 2J6 + ... */
/* is formed in . Use is then made of the property that */
/* 1 = J0 + 2J2 + 2J4 + 2J6 + ... to scale the resulting Js, */
/* by dividing the Js by . */
/* As the trial values become very large, scaling is applied */
/* when a J value held in rec2 exceeds . */
Rec: rec0 = 0; rec1 = 1; sum = 0;
Threshold = HUGE(rec2)/1.00000E5 * ABS(x/k);
DO p = k TO 1 BY -1;
rec2 = 2*p/x * rec1 - rec0;
IF p > n+1 THEN
J(n) = rec2;
ELSE
J(p-1) = rec2;
IF p = 1 THEN
sum = sum + rec2;
ELSE IF MOD(p, 2) = 1 THEN
sum = sum + 2*rec2;
IF ABS(rec2) > Threshold THEN
DO; /* Scale the results so as to avoid huge values. */
rec1 = rec1/Threshold; rec2 = rec2/Threshold;
sum = sum/Threshold;
DO m = n TO p-1 BY -1;
J(m) = J(m)/Threshold;
END;
put skip list ('**eps=', eps);
END;
rec0 = rec1; rec1 = rec2;
END; /* of the application of the recurrence formula. */
/* Stage 2: Normalize the J results. */
DO p = 0 TO n;
J(p) = J(p) / sum;
END;
IF s THEN
DO;
max = 0;
DO p = 0 TO n;
err = ABS(J(p) - Jbar(p));
IF err > max THEN max = err;
END; /* maximum error. */
IF max <= eps THEN RETURN;
put skip list ('EPS=', eps, ' max=', max, ' k=', k);
END;
ELSE
s = '1'b;
DO p = 0 TO n;
Jbar(p) = J(p);
END;
k = k + dist;
GO TO Rec;
END BESSELSETINT;