/************************************************************************/ /* */ /* 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;