ccccccccccccccccccccccccc     Program 5.5     cccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE BSSL (BJ,BY,N,X)
C
C Subroutine to generate J_n(x) and Y_n(x) with given x and
C up to N=NMAX-NTEL.
C 
      PARAMETER(NMAX = 30)
      PARAMETER(NTEL =  5)
      DIMENSION BJ(0:N),BY(0:N),B1(0:NMAX)
C
      IF(N.GT.(NMAX-NTEL)) STOP 'N is too large'
C
      PI = 4.0*ATAN(1.)
      GAMMA = 0.5772156649
C
      B1(NMAX)   = 0.0
      B1(NMAX-1) = 1.0
C
      SUM = 0.0
      DO   100  I  = NMAX-1, 1, -1 
        B1(I-1) = 2*I*B1(I)/X-B1(I+1)
        IF(MOD(I,2).EQ.0) SUM = SUM+2.0*B1(I)
  100 CONTINUE
      SUM = SUM + B1(0)
C
      DO   200  I  = 0, N
        BJ(I) = B1(I)/SUM
  200 CONTINUE
C
C Generating Y_n(x) starts here
C
      SUM1 = 0.0
      DO   300  K  = 1, NMAX/2
        SUM1 =SUM1+(-1)**K*B1(2*K)/K
  300 CONTINUE
C
      SUM1 = -4.0*SUM1/(PI*SUM)
      BY(0) = 2.0*(ALOG(X/2.0)+GAMMA)*BJ(0)/PI+SUM1
      BY(1) = (BJ(1)*BY(0)-2.0/(PI*X))/BJ(0)
C
      DO   400  I  = 1, N-1
        BY(I+1) = 2*I*BY(I)/X-BY(I-1)
  400 CONTINUE
C
      RETURN
      END