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