SUBROUTINE B1SLR(X,NB,IZE,B,NCALC) C THIS ROUTINE CALCULATES BESSEL FUNCTIONS I AND J OF REAL C ARGUMENT AND INTEGER ORDER. C C C EXPLANATION OF VARIABLES IN THE CALLING SEQUENCE C C X REAL ARGUMENT FOR WHICH I*S OR J*S C ARE TO BE CALCULATED. IF I*S ARE TO BE CALCULATED, C ABS(X) MUST BE LESS THAN EXPARG (WHICH SEE BELOW). C NB INTEGER TYPE. 1 + HIGHEST ORDER TO BE CALCULATED. C IT MUST BE POSITIVE. C IZE INTEGER TYPE. ZERO IF J*S ARE TO BE CALCULATED, 1 C IF I*S ARE TO BE CALCULATED. C B REAL VECTOR OF LENGTH NB, NEED NOT BE C INITIALIZED BY USER. IF THE ROUTINE TERMINATES C NORMALLY (NCALC=NB), IT RETURNS J(OR I)-SUB-ZERO C THROUGH J(OR I)-SUB-NB-MINUS-ONE OF X IN THIS C VECTOR. C NCALC INTEGER TYPE, NEED NOT BE INITIALIZED BY USER. C BEFORE USING THE RESULTS, THE USER SHOULD CHECK THAT C NCALC=NB, I.E. ALL ORDERS HAVE BEEN CALCULATED TO C THE DESIRED ACCURACY. SEE ERROR RETURNS BELOW. C C C EXPLANATION OF MACHINE-DEPENDENT CONSTANTS C C NSIG DECIMAL SIGNIFICANCE DESIRED. SHOULD BE SET TO C IFIX(ALOG10(2)*NBIT+1), WHERE NBIT IS THE NUMBER OF C BITS IN THE MANTISSA OF A REAL VARIABLE. C SETTING NSIG LOWER WILL RESULT IN DECREASED ACCURACY C WHILE SETTING NSIG HIGHER WILL INCREASE CPU TIME C WITHOUT INCREASING ACCURACY. THE TRUNCATION ERROR C IS LIMITED TO T=.5*10**-NSIG FOR J*S OF ORDER LESS C THAN ARGUMENT, AND TO A RELATIVE ERROR OF T FOR C I*S AND THE OTHER J*S. C NTEN LARGEST INTEGER K SUCH THAT 10**K IS MACHINE- C REPRESENTABLE IN SINGLE PRECISION. C LARGEX UPPER LIMIT ON THE MAGNITUDE OF X. BEAR IN MIND C THAT IF ABS(X)=N, THEN AT LEAST N ITERATIONS OF THE C BACKWARD RECURSION WILL BE EXECUTED. C EXPARG LARGEST REAL ARGUMENT THAT THE LIBRARY C EXP ROUTINE CAN HANDLE. C C PORT NOTE, SEPTEMBER 8,1976 - C THE LARGEX AND EXPARG TESTS ARE MADE IN THE OUTER ROUTINES - C BESRJ AND BESRI, WHICH CALL B1SLR. C C C ERROR RETURNS C C PORT NOTE, SEPTEMBER 8, 1976 - C THE NOTES BELOW ARE KEPT IN FOR THE RECORD, BUT, AS ABOVE, C THE ACTUAL TESTS ARE NOW IN THE OUTER CALLING ROUTINES. C C LET G DENOTE EITHER I OR J. C IN CASE OF AN ERROR, NCALC.NE.NB, AND NOT ALL G*S C ARE CALCULATED TO THE DESIRED ACCURACY. C IF NCALC.LT.0, AN ARGUMENT IS OUT OF RANGE. NB.LE.0 C OR IZE IS NEITHER 0 NOR 1 OR IZE=1 AND ABS(X).GE.EXPARG. C IN THIS CASE, THE B-VECTOR IS NOT CALCULATED, AND NCALC C IS SET TO MIN0(NB,0)-1 SO NCALC.NE.NB. C NB.GT.NCALC.GT.0 WILL OCCUR IF NB.GT.MAGX AND ABS(G- C SUB-NB-OF-X/G-SUB-MAGX+NP-OF-X).LT.10.**(NTEN/2), I.E. NB C IS MUCH GREATER THAN MAGX. IN THIS CASE, B(N) IS CALCU- C LATED TO THE DESIRED ACCURACY FOR N.LE.NCALC, BUT FOR C NCALC.LT.N.LE.NB, PRECISION IS LOST. IF N.GT.NCALC AND C ABS(B(NCALC)/B(N)).EQ.10**-K, THEN THE LAST K SIGNIFICANT C FIGURES OF B(N) ARE ERRONEOUS. IF THE USER WISHES TO C CALCULATE B(N) TO HIGHER ACCURACY, HE SHOULD USE AN C ASYMPTOTIC FORMULA FOR LARGE ORDER. C REAL 1 X,B,P,TEST,TEMPA,TEMPB,TEMPC,SIGN,SUM,TOVER, 2 PLAST,POLD,PSAVE,PSAVEL,R1MACH DIMENSION B(NB) DATA NSIG/0/, NTEN/0/ IF(NSIG .NE. 0) GO TO 1 NSIG = IFIX(-ALOG10(R1MACH(3))+1.) NTEN = ALOG10(R1MACH(2)) 1 TEMPA=ABS(X) MAGX=IFIX((TEMPA)) C SIGN=FLOAT(1-2*IZE) NCALC=NB C USE 2-TERM ASCENDING SERIES FOR SMALL X IF(TEMPA**4.LT..1E0**NSIG) GO TO 30 C INITIALIZE THE CALCULATION OF P*S NBMX=NB-MAGX N=MAGX+1 PLAST=1.E0 P=FLOAT(2*N)/TEMPA C CALCULATE GENERAL SIGNIFICANCE TEST TEST=2.E0*1.E1**NSIG IF(IZE.EQ.1.AND.2*MAGX.GT.5*NSIG) TEST=SQRT(TEST*P) IF(IZE.EQ.1.AND.2*MAGX.LE.5*NSIG) TEST=TEST/1.585**MAGX M=0 IF(NBMX.LT.3) GO TO 4 C CALCULATE P*S UNTIL N=NB-1. CHECK FOR POSSIBLE OVERFLOW. TOVER=1.E1**(NTEN-NSIG) NSTART=MAGX+2 NEND=NB-1 DO 3 N=NSTART,NEND POLD=PLAST PLAST=P P=FLOAT(2*N)*PLAST/TEMPA-SIGN*POLD IF(P-TOVER) 3,3,5 3 CONTINUE C CALCULATE SPECIAL SIGNIFICANCE TEST FOR NBMX.GT.2. TEST=AMAX1(TEST,SQRT(PLAST*1.E1**NSIG)*SQRT(2.E0*P)) C CALCULATE P*S UNTIL SIGNIFICANCE TEST PASSES 4 N=N+1 POLD=PLAST PLAST=P P=FLOAT(2*N)*PLAST/TEMPA-SIGN*POLD IF(P.LT.TEST) GO TO 4 IF(IZE.EQ.1.OR.M.EQ.1) GO TO 12 C FOR J*S, A STRONG VARIANT OF THE TEST IS NECESSARY. C CALCULATE IT, AND CALCULATE P*S UNTIL THIS TEST IS PASSED. M=1 TEMPB=P/PLAST TEMPC=FLOAT(N+1)/TEMPA IF(TEMPB+1.E0/TEMPB.GT.2.E0*TEMPC)TEMPB=TEMPC+SQRT(TEMPC**2-1.E0) TEST=TEST/SQRT(TEMPB-1.E0/TEMPB) IF(P-TEST) 4,12,12 C TO AVOID OVERFLOW, DIVIDE P*S BY TOVER. CALCULATE P*S C UNTIL ABS(P).GT.1. 5 TOVER=1.E1**NTEN P=P/TOVER PLAST=PLAST/TOVER PSAVE=P PSAVEL=PLAST NSTART=N+1 6 N=N+1 POLD=PLAST PLAST=P P=FLOAT(2*N)*PLAST/TEMPA-SIGN*POLD IF(P.LE.1.E0) GO TO 6 TEMPB=FLOAT(2*N)/TEMPA IF(IZE.EQ.1) GO TO 8 TEMPC=.5E0*TEMPB TEMPB=PLAST/POLD IF(TEMPB+1.E0/TEMPB.GT.2.E0*TEMPC)TEMPB=TEMPC+SQRT(TEMPC**2-1.E0) C CALCULATE BACKWARD TEST, AND FIND NCALC, THE HIGHEST N C SUCH THAT THE TEST IS PASSED. 8 TEST=.5E0*POLD*PLAST*(1.E0-1.E0/TEMPB**2)/1.E1**NSIG P=PLAST*TOVER N=N-1 NEND=MIN0(NB,N) DO 9 NCALC=NSTART,NEND POLD=PSAVEL PSAVEL=PSAVE PSAVE=FLOAT(2*N)*PSAVEL/TEMPA-SIGN*POLD IF(PSAVE*PSAVEL-TEST) 9,9,10 9 CONTINUE NCALC=NEND+1 10 NCALC=NCALC-1 C THE SUM B(1)+2B(3)+2B(5)... IS USED TO NORMALIZE. M, THE C COEFFICIENT OF B(N), IS INITIALIZED TO 2 OR 0. 12 N=N+1 M=2*N-4*(N/2) C INITIALIZE THE BACKWARD RECURSION AND THE NORMALIZATION C SUM TEMPB=0.E0 TEMPA=1.E0/P SUM=FLOAT(M)*TEMPA NEND=N-NB IF(NEND) 17,15,13 C RECUR BACKWARD VIA DIFFERENCE EQUATION, CALCULATING (BUT C NOT STORING) B(N), UNTIL N=NB. 13 DO 14 L=1,NEND N=N-1 TEMPC=TEMPB TEMPB=TEMPA TEMPA=FLOAT(2*N)*TEMPB/X-SIGN*TEMPC M=2-M 14 SUM=SUM+FLOAT(M)*TEMPA C STORE B(NB) 15 B(N)=TEMPA IF(NB.GT.1) GO TO 16 C NB=1. SINCE 2*TEMPA WAS ADDED TO THE SUM, TEMPA MUST BE C SUBTRACTED SUM=SUM-TEMPA GO TO 23 C CALCULATE AND STORE B(NB-1) 16 N=N-1 B(N) =FLOAT(2*N)*TEMPA/X-SIGN*TEMPB IF(N.EQ.1) GO TO 22 M=2-M SUM=SUM+FLOAT(M)*B(N) GO TO 19 C N.LT.NB, SO STORE B(N) AND SET HIGHER ORDERS TO ZERO 17 B(N)=TEMPA NEND=-NEND DO 18 L=1,NEND K=N+L 18 B(K)=0.E0 19 NEND=N-2 IF(NEND.EQ.0) GO TO 21 C CALCULATE VIA DIFFERENCE EQUATION AND STORE B(N), C UNTIL N=2 DO 20 L=1,NEND N=N-1 B(N)=(FLOAT(2*N)*B(N+1))/X-SIGN*B(N+2) M=2-M 20 SUM=SUM+FLOAT(M)*B(N) C CALCULATE B(1) 21 B(1)=2.E0*B(2)/X-SIGN*B(3) 22 SUM=SUM+B(1) C NORMALIZE--IF IZE=1, DIVIDE SUM BY COSH(X). DIVIDE ALL C B(N) BY SUM. 23 IF(IZE.EQ.0) GO TO 25 TEMPA=EXP(ABS(X)) SUM=2.E0*SUM/(TEMPA+1.E0/TEMPA) 25 DO 26 N=1,NB 26 B(N)=B(N)/SUM RETURN C C TWO-TERM ASCENDING SERIES FOR SMALL X 30 TEMPA=1.E0 TEMPB=-.25E0*X*X*SIGN B(1)=1.E0+TEMPB IF(NB.EQ.1) GO TO 32 DO 31 N=2,NB TEMPA=TEMPA*X/FLOAT(2*N-2) 31 B(N)=TEMPA*(1.E0+TEMPB/FLOAT(N)) 32 RETURN END