ADAS Subroutine rd2bs

       FUNCTION RD2BS(N,L,N2,L2)
       IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C
C  PURPOSE: GENERATION OF HYDROGENIC BOUND-BOUND RADIAL INTEGRALS USING
C           RECURRENCE RELATIONS.
C
C-----------------------------------------------------------------------
C VERSION  : 1.1
C DATE     : 18-03-1999
C MODIFIED : ???
C
C VERSION  : 1.2
C DATE     : 05-10-2000
C MODIFIED : ???
C             - Removed junk from columns > 72
C
C VERSION  : 1.3
C DATE     : 16-05-2007
C MODIFIED : Allan Whiteford
C             - Updated comments as part of subroutine documentation
C               procedure.
C
C-----------------------------------------------------------------------
       SC=64.0
       SCL=0.015625
       EN=N
       EN2=EN*EN
       EK=N2
       EK=1.0/EK
       EK2=-EK*EK
       V=1.0+EN2*EK2
       W=1.0+EN*EK
       W=W*W
       U=8.0*EN2/(W*W)
       P=1.0
       JS=0
       SC2=SC*SC
       SCL2=SCL*SCL
       DO 5 I=1,N
       EI=I
       P=P*U*(1.0+EI*EI*EK2)/(EI*(2.0*EI-1.0))
       AP=DABS(P)
       IF(SCL2.LE.AP)GO TO 5
       JS=JS-1
       P=SC2*P
    5  CONTINUE
       P=2.0*EN*EK*P
       P=4.0*EK*DSQRT(P)/(V*V)
       NL=N+1
       DO 10 I=NL,N2
       P=P*V/W
       AP=DABS(P)
       IF(SCL.LE.AP)GO TO 10
       JS=JS-1
       P=SC*P
   10  CONTINUE
       IF(L2.EQ.L+1)GO TO 11
       IF(L2.EQ.L-1)GO TO 20
       RD2BS=0.0
       GO TO 7
   11  T2=P
       U=(2.0*EN-1.0)*V
       U=DSQRT(U)
       T3=0.5*U*T2
       NU=N-2
       IF(L-NU)14,13,12
   12  T3=T2
   13  GO TO 40
   14  DO 16 I=L2,NU
       LI=NU-I+L
       EL1=LI+1
       EL2=LI+2
       ES=EL2*EL2
       T1=T2
       T2=T3
       T3=(4.0*(EN2-ES)+EL2*(2.0*EL2-1.0)*V)*T2-2.0*EN*U*T1
       U=(EN2-EL1*EL1)*(1.0+ES*EK2)
       U=DSQRT(U)
       T3=T3/(2.0*EN*U)
       AT3=DABS(T3)
       IF(AT3.LE.SC)GO TO 16
       JS=JS+1
       T3=SCL*T3
       T2=SCL*T2
   16  CONTINUE
       GO TO 40
   20  T2=P
       EN1=N-1
       U=V/(1.0+EN1*EN1*EK2)
       T2=DSQRT(U)*T2/(2.0*EN)
       U=(2.0*EN-1.0)*(1.0+(EN-2.0)*(EN-2.0)*EK2)
       U=DSQRT(U)
       T3=(4.0+EN1*V)*(2.0*EN-1.0)*T2/(2.0*EN*U)
       NU=N-3
       IF(L-NU-1)24,23,22
   22  T3=T2
   23  GO TO 40
   24  DO 26 I=L,NU
       LI=NU-I+L
       EL=LI
       EL1=LI+1
       ES=EL1*EL1
       T1=T2
       T2=T3
       T3=(4.0*(EN2-ES)+EL1*(2.0*EL1+1.0)*V)*T2-2.0*EN*U*T1
       U=(EN2-ES)*(1.0+EL*EL*EK2)
       U=DSQRT(U)
       T3=T3/(2.0*EN*U)
       AT3=DABS(T3)
       IF(AT3.LE.SC)GO TO 26
       JS=JS+1
       T3=SCL*T3
       T2=SCL*T2
   26  CONTINUE
   40  RD2BS=EN2*EN2*T3*T3*4096.0**JS
    7  RETURN
      END
      INTEGER             L,           L2,          N,           N2