(in-package :f77)

(def-fortran-routine dbesj :void
"***BEGIN PROLOGUE  DBESJ
***PURPOSE  Compute an N member sequence of J Bessel functions
            J/SUB(ALPHA+K-1)/(X), K=1,...,N for non-negative ALPHA
            and X.
***LIBRARY   SLATEC
***CATEGORY  C10A3
***TYPE      DOUBLE PRECISION (BESJ-S, DBESJ-D)
***KEYWORDS  J BESSEL FUNCTION, SPECIAL FUNCTIONS
***AUTHOR  Amos, D. E., (SNLA)
           Daniel, S. L., (SNLA)
           Weston, M. K., (SNLA)
***DESCRIPTION

     Abstract  **** a double precision routine ****
         DBESJ computes an N member sequence of J Bessel functions
         J/sub(ALPHA+K-1)/(X), K=1,...,N for non-negative ALPHA and X.
         A combination of the power series, the asymptotic expansion
         for X to infinity and the uniform asymptotic expansion for
         NU to infinity are applied over subdivisions of the (NU,X)
         plane.  For values of (NU,X) not covered by one of these
         formulae, the order is incremented or decremented by integer
         values into a region where one of the formulae apply. Backward
         recursion is applied to reduce orders by integer values except
         where the entire sequence lies in the oscillatory region.  In
         this case forward recursion is stable and values from the
         asymptotic expansion for X to infinity start the recursion
         when it is efficient to do so. Leading terms of the series and
         uniform expansion are tested for underflow.  If a sequence is
         requested and the last member would underflow, the result is
         set to zero and the next lower order tried, etc., until a
         member comes on scale or all members are set to zero.
         Overflow cannot occur.

         The maximum number of significant digits obtainable
         is the smaller of 14 and the number of digits carried in
         double precision arithmetic.

     Description of Arguments

         Input      X,ALPHA are double precision
           X      - X .GE. 0.0D0
           ALPHA  - order of first member of the sequence,
                    ALPHA .GE. 0.0D0
           N      - number of members in the sequence, N .GE. 1

         Output     Y is double precision
           Y      - a vector whose first N components contain
                    values for J/sub(ALPHA+K-1)/(X), K=1,...,N
           NZ     - number of components of Y set to zero due to
                    underflow,
                    NZ=0   , normal return, computation completed
                    NZ .NE. 0, last NZ components of Y set to zero,
                             Y(K)=0.0D0, K=N-NZ+1,...,N.

     Error Conditions
         Improper input arguments - a fatal error
         Underflow  - a non-fatal error (NZ .NE. 0)

***REFERENCES  D. E. Amos, S. L. Daniel and M. K. Weston, CDC 6600
                 subroutines IBESS and JBESS for Bessel functions
                 I(NU,X) and J(NU,X), X .GE. 0, NU .GE. 0, ACM
                 Transactions on Mathematical Software 3, (1977),
                 pp. 76-92.
               F. W. J. Olver, Tables of Bessel Functions of Moderate
                 or Large Orders, NPL Mathematical Tables 6, Her
                 Majesty's Stationery Office, London, 1962.
"
  (X :double-float :input)
  (alpha :double-float :input)
  (n :integer :input)
  (Y (* :double-float) :output)
  (NZ :integer :output))