SUBROUTINE DB1SLR(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     DOUBLE PRECISION 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     DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE 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 DOUBLE PRECISION ARGUMENT THAT THE LIBRARY
C       DEXP ROUTINE CAN HANDLE.
C
C PORT NOTE, SEPTEMBER 8,1976 -
C THE LARGEX AND EXPARG TESTS ARE MADE IN THE OUTER ROUTINES -
C DBESRJ AND DBESRI, WHICH CALL DB1SLR.
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
      DOUBLE PRECISION DLOG10, DSQRT, DEXP,
     1 X,B,P,TEST,TEMPA,TEMPB,TEMPC,SIGN,SUM,TOVER,
     2 PLAST,POLD,PSAVE,PSAVEL,D1MACH
      DIMENSION B(NB)
      DATA NSIG/0/, NTEN/0/
      IF(NSIG .NE. 0) GO TO 1
      NSIG = IFIX(-ALOG10(SNGL(D1MACH(3)))+1.)
      NTEN = DLOG10(D1MACH(2))
    1 TEMPA=DABS(X)
      MAGX=IFIX(SNGL(TEMPA))
C
      SIGN=DBLE(FLOAT(1-2*IZE))
      NCALC=NB
C USE 2-TERM ASCENDING SERIES FOR SMALL X
      IF(TEMPA**4.LT..1D0**NSIG) GO TO 30
C INITIALIZE THE CALCULATION OF P*S
      NBMX=NB-MAGX
      N=MAGX+1
      PLAST=1.D0
      P=DBLE(FLOAT(2*N))/TEMPA
C CALCULATE GENERAL SIGNIFICANCE TEST
      TEST=2.D0*1.D1**NSIG
      IF(IZE.EQ.1.AND.2*MAGX.GT.5*NSIG) TEST=DSQRT(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.D1**(NTEN-NSIG)
      NSTART=MAGX+2
      NEND=NB-1
      DO 3 N=NSTART,NEND
      POLD=PLAST
      PLAST=P
      P=DBLE(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=DMAX1(TEST,DSQRT(PLAST*1.D1**NSIG)*DSQRT(2.D0*P))
C CALCULATE P*S UNTIL SIGNIFICANCE TEST PASSES
    4 N=N+1
      POLD=PLAST
      PLAST=P
      P=DBLE(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=DBLE(FLOAT(N+1))/TEMPA
      IF(TEMPB+1.D0/TEMPB.GT.2.D0*TEMPC)TEMPB=TEMPC+DSQRT(TEMPC**2-1.D0)
      TEST=TEST/DSQRT(TEMPB-1.D0/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.D1**NTEN
      P=P/TOVER
      PLAST=PLAST/TOVER
      PSAVE=P
      PSAVEL=PLAST
      NSTART=N+1
    6 N=N+1
      POLD=PLAST
      PLAST=P
      P=DBLE(FLOAT(2*N))*PLAST/TEMPA-SIGN*POLD
      IF(P.LE.1.D0) GO TO 6
      TEMPB=DBLE(FLOAT(2*N))/TEMPA
      IF(IZE.EQ.1) GO TO 8
      TEMPC=.5D0*TEMPB
      TEMPB=PLAST/POLD
      IF(TEMPB+1.D0/TEMPB.GT.2.D0*TEMPC)TEMPB=TEMPC+DSQRT(TEMPC**2-1.D0)
C CALCULATE BACKWARD TEST, AND FIND NCALC, THE HIGHEST N
C SUCH THAT THE TEST IS PASSED.
    8 TEST=.5D0*POLD*PLAST*(1.D0-1.D0/TEMPB**2)/1.D1**NSIG
      P=PLAST*TOVER
      N=N-1
      NEND=MIN0(NB,N)
      DO 9 NCALC=NSTART,NEND
      POLD=PSAVEL
      PSAVEL=PSAVE
      PSAVE=DBLE(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.D0
      TEMPA=1.D0/P
      SUM=DBLE(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=(DBLE(FLOAT(2*N))*TEMPB)/X-SIGN*TEMPC
      M=2-M
   14 SUM=SUM+DBLE(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) =(DBLE(FLOAT(2*N))*TEMPA)/X-SIGN*TEMPB
      IF(N.EQ.1) GO TO 22
      M=2-M
      SUM=SUM+DBLE(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.D0
   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)=(DBLE(FLOAT(2*N))*B(N+1))/X-SIGN*B(N+2)
      M=2-M
   20 SUM=SUM+DBLE(FLOAT(M))*B(N)
C CALCULATE B(1)
   21 B(1)=2.D0*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=DEXP(DABS(X))
      SUM=2.D0*SUM/(TEMPA+1.D0/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.D0
      TEMPB=-.25D0*X*X*SIGN
      B(1)=1.D0+TEMPB
      IF(NB.EQ.1) GO TO 32
      DO 31 N=2,NB
      TEMPA=TEMPA*X/DBLE(FLOAT(2*N-2))
   31 B(N)=TEMPA*(1.D0+TEMPB/DBLE(FLOAT(N)))
   32 RETURN
      END