      subroutine DBESJN(X,ALPHA,NUM,BJ)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c>> 1996-03-30 DBESJN Krogh  Removed INT from type statement.
C>> 1995-11-13 DBESJN Krogh  Converted SFTRAN to Fortran
C>> 1995-10-24 DBESJN Krogh  Removed blanks in numbers for C conversion.
C>> 1994-10-19 DBESJN Krogh  Changes to use M77CON
C>> 1994-04-19 DBESJN CLL  Edited to make DP & SP files similar.
C>> 1992-03-13 DBESJN FTK  Removed implicit statements.
C>> 1989-08-09 DBESJN CLL  More accurate HALFPI for Cray
C>> 1986-03-18 DBESJN Lawson  Initial code.
c--D replaces "?": ?BESJN, ?GAMMA, ?BESPQ, ?ERV1
C
C     COMPUTE THE J-BESSEL FUNCTIONS OF X FOR ORDER ALPHA THROUGH
C     ALPHA + NUM - 1 in steps of one.
C     STORE THE RESULT INTO BJ(I),I=1,...,NUM.
C     Require X .ge. 0, ALPHA .ge. 0, NUM .ge. 1
C
c     Original code due to E. W. Ng and W. V. Snyder, JPL, 1973.
c     Modified by Ng and S. Singletary, 1974.
C     C.Lawson & S.Chan, JPL, 1984 Apr 05:
c        Changed calling sequence.
c        Adapted to SFTRAN3 and Fortran 77.
c        Added code to avoid overflow during the recurrsion.
c        Added code to use Taylor series for small x.
c        March 23. Improved accuracy of backward recursion.
c        Changed to use P and Q instead of R and THETA
c                in the region of large X.
C     ------------------------------------------------------------------
C
      save EPS, SMALL, XPQ, BIG, HICUT
c     ----------
      external           DGAMMA, DBESPQ, D1MACH, ERMSG, DERV1, IERV1
c     ----------
      integer I,I1,II,J,K,M,MU,NMAX,NMIN,NUM
      double precision D1MACH, DGAMMA
      double precision ALPHA,BESV,BESVM1,BESVP1,BIG,BJ(NUM),BJNU,BOT
      double precision C11293,C16,C1P5,C59,CHI,CP6
      double precision DR,EM,EMU,EPS,ETA,FAC,FAC2,FK,FKM1,FKP1,FV,FVM1
      double precision FVP1,G,GNU,HALF,HALFPI,HALFX,HICUT
      double precision ONE,P,Q,SCALE,SMALL,SUM,T1,T2,TENTH,TERM
      double precision TERM1,TEST,TOP,TWO,TWODX,V,VPMU,X,XPQ,ZERO
      parameter(ZERO=0.D0, ONE=1.D0, TWO=2.D0, C16=16.D0)
      parameter( TENTH = 0.1D0, HALF = 0.5D0)
      parameter(HALFPI = 1.57079632679489661923132169163975144D0)
      parameter(C11293 = 1.1293D0, CP6 = 0.60206D0, C59 = 0.59D0)
      parameter(C1P5 = 1.5D0)
      data EPS / ZERO /
C     ------------------------------------------------------------------
C
c     Set environment parameters.  This should happen only
c     the first time this subroutine is called.
c
      if( EPS .EQ. ZERO ) then
         EPS = D1MACH(3)
         HICUT = ONE /(EPS*C16)
         XPQ = C11293 * (CP6 - log10(EPS)) - C59
         SMALL = C16 * D1MACH(1) / EPS
         BIG = D1MACH(2)/TWO
      end if
C
c     Test validity of input values.
c
      if (X .lt. ZERO .or. ALPHA .lt. ZERO .or. NUM .lt. 1) then
c
c                                                           Error 1.
        call ERMSG('DBESJN',1,0,
     *             'REQUIRE X.GE.0, ALPHA.GE.0, NUM.GE.1',',')
      else
c                            Begin computation.
         NMIN = int(ALPHA)
         V = ALPHA - dble(NMIN)
         NMAX = NMIN + NUM - 1
C
         if( X .le. TENTH ) then
c
c ********************* Code for the small X case. *********************
c
            if( X .eq. ZERO ) then
c                                        Special for X .eq. 0.
               do 80 I = 1, NUM
                  BJ(I) =  ZERO
   80          continue
               if (ALPHA .eq. 0) BJ(1) = ONE
            else
c                    Here use Taylor series for small x.
C
               GNU = ALPHA
               HALFX = HALF*X
               FAC2 = -HALFX*HALFX
               TERM1 = (HALFX**GNU) / DGAMMA(GNU + ONE)
               do 150 I = 1, NUM
c           Sum the series for the Bessel fcn J sub GNU of X given
c           in Eq 9.1.10, page 360, of AMS 55.
c           1984 March 9, JPL, C. L. Lawson.
c
                  if(TERM1 .eq. ZERO) then
                     BJNU = ZERO
                  else
c
                     SUM = ZERO
                     TOP = FAC2
                     BOT = GNU + ONE
                     T1 = ONE
                     T2 = BOT
                     TERM = TOP/BOT
c
  100                if ( abs(TERM) .GT. EPS ) then
                        SUM = SUM + TERM
                        TOP = TOP * FAC2
                        T1 = T1 + ONE
                        T2 = T2 + ONE
                        BOT = BOT * T1 * T2
                        TERM = TOP/BOT
                        go to 100
                     end if
                     BJNU = TERM1 + TERM1 * SUM
                  end if
c
                  BJ(I) = BJNU
                  if( BJNU .eq. ZERO ) then
C
c     Here current result has underflowed to zero, so we will
c     set the rest of the results to zero also.
c
                     do 120 J = I+1, NUM
                        BJ(J) = ZERO
  120                continue
                     return
                  end if
                  GNU = GNU + ONE
                  TERM1 = TERM1 * (HALFX / GNU)
  150          continue
            end if
            return
         else if( X .lt. MAX(dble(NMAX+1), XPQ) ) then
c
c ********************* Code for the middle X case. ********************
c
C
C     J-TYPE BESSEL FUNCTIONS FOLLOW THE RECURRENCE RELATION
C     F(V-1,X)=(2*V/X)*F(V,X)-F(V+1,X).
C
            TWODX = TWO / X
            MU = MAX( int(X)+1, NMAX)
            DR = TWODX * (V+dble(MU))
            FKP1 = ONE
            FK =  ZERO
C
C     RECUR FORWARD UNTIL FKP1 IS GREATER THAN PRECISION OF ARITHMETIC.
C
  200       if ( EPS * abs(FKP1) .le. ONE ) then
               MU = MU + 1
               DR = DR + TWODX
               FKM1 = FK
               FK = FKP1
               FKP1 = DR * FK - FKM1
               go to 200
            end if
C
C     WE ARE NOW ASSURED THAT BACKWARD RECURRENCE FROM MU WILL YIELD
C     ACCURATE RESULTS.
C
C                                        GUARANTEE EVEN MU
            if (MOD(MU,2) .ne. 0)  MU = MU + 1
            FVM1 = SMALL
            FV = ZERO
            ETA = ONE
            SUM = FVM1
            M = MU / 2
            EM = dble(M)
            EMU = dble(MU)
            FAC = (V + EMU) * TWODX
C
c     Set TEST = largest value that can be multiplied by
c     FAC without risking overflow.  The present value of
c     FAC is the largest that will occur during the recursion.
c     TEST will be used to protect against overflow during
c     the recursion.
c
            TEST = BIG / MAX(ONE, FAC)
C
C                            Loop while MU .GT. ZERO
  220       continue
               FVP1 = FV
               FV = FVM1
               if( abs(FV) .GT. TEST )then
                  FV = FV / SUM
                  FVP1 = FVP1 / SUM
                  I1 = MAX( 1, MU - NMIN + 1 )
                  do 230 II = I1, NUM
                     BJ(II) = BJ(II) / SUM
  230             continue
                  SUM = ONE
               end if
               FVM1 = FAC * FV - FVP1
               MU = MU -1
               EMU = EMU - ONE
               FAC = (V + EMU) * TWODX
               if (MU .ge. NMIN .AND. MU .le. NMAX) BJ(MU-NMIN+1) = FVM1
               if (MOD(MU,2) .eq. 0) then
                  if (V .eq. ZERO)  then
                     SUM = SUM + FVM1
                     if (MU .eq. 0) then
                        SCALE = ONE / SUM
                        go to 250
                     end if
                     SUM = SUM + FVM1
                  else
                     if (MU .ne. 0) then
                        VPMU = V + EMU
                        ETA = ETA * (EM/(V+(EM-ONE)))*(VPMU/(VPMU+TWO))
                        SUM = SUM + FVM1 * ETA
                        EM = EM - ONE
                     else
c
c           Here MU = 0 and EM = 0NE.  Thus the expression for
c           updating ETA reduces to the following simpler
c           expression.
c
                        ETA = ETA / (V + TWO)
                        SUM = SUM + FVM1 * ETA
                        SCALE = DGAMMA(V+ONE) / ETA * SUM * TWODX ** V
                        SCALE = ONE / SCALE
                        go to 250
                     end if
                  end if
               end if
            go to 220
  250       continue
C
C     NORMALIZE BJ() TO GET VALUES OF J-BESSEL FUNCTION.
C
            do 260 I = 1, NUM
               BJ(I) = BJ(I) * SCALE
  260       continue
            return
         else if( X .lt. HICUT ) then
c
c ********************* Code for the large X case. *********************
c
c  Here we have X .ge. XPQ, and V in [0.,1.).
c     The asymptotic series for  the auxiliary functions P and Q can be
c     used.  From these we will compute J(V,X) and J(V+1,X) and
c     then recur forward.  Reference: NBS AMS 55 Eqs 9.2.5 & 9.2.6
c
            call DBESPQ (X,V,  P,Q)
            CHI = X - (V + HALF) * HALFPI
            BESV = sqrt(ONE / (HALFPI*X)) * (P*COS(CHI) - Q*SIN(CHI))
C
            if( NMAX .GT. 0 ) then
               call DBESPQ (X,V+ONE,   P,Q)
               CHI = X - (V + C1P5) * HALFPI
               BESVP1 = sqrt(ONE / (HALFPI*X))*(P*COS(CHI)-Q*SIN(CHI))
            end if
c
            TWODX = TWO / X
c
c           Given BESV = J(V,X), BESVP1 = J(V+1,X), TWODX = 2/X,
c           NMIN, NUM, NMAX = NMIN + NUM -1, X, ALPHA, and BIG.
c           Recur forward and store J(NMIN+V) thru J(NMAX+V) in
c           BJ(1) thru BJ(NUM).
c           There should be no overflow posibility in this forward
c           recursion since NMAX .le. X - 1, and in this region the
c           magnitude of the J function is less than one.
c
            if( NMIN .eq. 0 ) then
               BJ(1) = BESV
               if( NMAX .GT. 0 ) then
                  BJ(2) = BESVP1
               end if
            else if( NMIN .eq. 1 ) then
               BJ(1) = BESVP1
            end if
c
            if( NMAX .GT. 1 ) then
               G = V * TWODX
c
c        Note:  In the following statement, 3-NMIN can be nonpositive.
c
               do 300 K = 3-NMIN, NUM
                  BESVM1 = BESV
                  BESV   = BESVP1
                  G = G + TWODX
                  BESVP1 = G * BESV - BESVM1
                  if( K .ge. 1)  BJ(K) = BESVP1
  300          continue
            end if
            return
         else
c                                                  Error 2.
            call ERMSG('DBESJN', 2, 0,
     *         'Cannot obtain any accuracy when X exceeds HICUT.', ',')
            call DERV1('HICUT', HICUT, ',')
         end if
      end if
      call DERV1('X',X,',')
      call DERV1('ALPHA',ALPHA,',')
      call IERV1('NUM',NUM,'.')
      return
      end
