subroutine ribesl ( x, alpha, nb, ize, b, ncalc )

!*****************************************************************************80
!
!! ribesl() computes bessel i functions with real (non integer) order.
!
!  Discussion:
!
!    this routine calculates bessel functions i sub(n+alpha) (x)
!    for non-negative argument x, and non-negative order n+alpha,
!    with or without exponential scaling.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
! explanation of variables in the calling sequence
!
! x     - working precision non-negative real argument for which
!         i's or exponentially scaled i's (i*exp(-x))
!         are to be calculated.  if i's are to be calculated,
!         x must be less than exparg (see below).
! alpha - working precision fractional part of order for which
!         i's or exponentially scaled i's (i*exp(-x)) are
!         to be calculated.  0 .le. alpha .lt. 1.0.
! nb    - integer number of functions to be calculated, nb .gt. 0.
!         the first function calculated is of order alpha, and the
!         last is of order (nb - 1 + alpha).
! ize   - integer type.  ize = 1 if unscaled i's are to calculated,
!         and 2 if exponentially scaled i's are to be calculated.
! b     - working precision output vector of length nb.  if the routine
!         terminates normally (ncalc=nb), the vector b contains the
!         functions i/alpha/(x) through i/nb-1+alpha/(x), or the
!         corresponding exponentially scaled functions.
! ncalc - integer output variable indicating possible errors.
!         before using the vector b, the user should check that
!         ncalc=nb, i.e., all orders have been calculated to
!         the desired accuracy.  see error returns below.
!
!
!*******************************************************************
!*******************************************************************
!
! explanation of machine-dependent constants
!
! nsig   - decimal significance desired.  should be set to
!          ifix(alog10(2)*nbit+1), where nbit is the number of
!          bits in the mantissa of a working precision variable.
!          setting nsig lower will result in decreased accuracy
!          while setting nsig higher will increase cpu time
!          without increasing accuracy.  the truncation error
!          is limited to a relative error of t=.5*10**(-nsig).
! enten  - 10.0 ** k, where k is the largest integer such that
!          enten is machine-representable in working precision.
! ensig  - 10.0 ** nsig.
! rtnsig - 10.0 ** (-k) for the smallest integer k such that
!          k .ge. nsig/4.
! enmten - the smallest abs(x) such that x/4 does not underflow.
! xlarge - upper limit on the magnitude of x when ize=2.  bear
!          in mind that if abs(x)=n, then at least n iterations
!          of the backward recursion will be executed.
! exparg - largest working precision argument that the library
!          exp routine can handle and upper limit on the
!          magnitude of x when ize=1.
!
!     approximate values for some important machines are:
!
!         ibm/195    cdc/7600  univac/1108   vax 11/780 (unix)
!          (d.p.)  (s.p.,rndg)    (d.p.)     (s.p.)     (d.p.)
!
! nsig      16         14          18           8        17
! enten   1.0d75     1.0e322     1.0d307     1.0e38    1.0d38
! ensig   1.0d16     1.0e14      1.0d18      1.0e8     1.0d17
! rtnsig  1.0d-4     1.0e-4      1.0d-5      1.0e-2    1.0d-4
! enmten  2.2d-78    1.0e-290    1.2d-308    1.2e-37   1.2d-37
! xlarge  1.0d4      1.0e4       1.0d4       1.0e4     1.0d4
! exparg  174.0d0    740.0e0     709.0d0     88.0e0    88.0d0
!
!*******************************************************************
!*******************************************************************
!
!
! error returns
!
!  in case of an error,  ncalc .ne. nb,  and not all i's are
!  calculated to the desired accuracy.
!
!  ncalc .lt. 0:  an argument is out of range. for example,
!     nb .le. 0, ize is not 1 or 2, or ize=1 and abs(x) .ge. exparg.
!     in this case, the b-vector is not calculated, and ncalc is
!     set to  min0(nb,0)-1  so that ncalc .ne. nb.
!
!  nb .gt. ncalc .gt. 0: not all requested function values could
!     be calculated accurately.  this usually occurs because nb is
!     much larger than abs(x).  in this case, b(n) is calculated
!     to the desired accuracy for  n .le. ncalc,  but precision
!     is lost for ncalc .lt. n .le. nb.  if b(n) does not vanish
!     for  n .gt. ncalc  (because it is too small to be represented),
!     and  b(n)/b(ncalc) = 10**(-k), then only the first nsig-k
!     significant figures of b(n) can be trusted.
!
!
! other subprograms required (single precision version)
!
!     exp,gamma,amax1,sqrt,float,ifix,min0
!
! other subprograms required (double precision version)
!
!     dble,dexp,dgamma,dmax1,dsqrt,float,ifix,min0,sngl
!
!
! acknowledgement
!
!  this program is based on a program written by david j.
!  sookne (2) that computes values of the bessel functions j or
!  i of real argument and integer order.  modifications include
!  the restriction of the computation to the i bessel function
!  of non-negative real argument, the extension of the computation
!  to arbitrary positive order, the inclusion of optional
!  exponential scaling, and the elimination of most underflow.
!
!  Author:
!
!    William Cody
!
!  Reference:
!
!  (1) olver, f. w. j., and sookne, d. j., "a note on backward
!        recurrence algorithms," math. comp. 26, 1972, pp 941-947.
!
!  (2) sookne, d. j., "bessel functions of real argument and
!        integer order," nbs jour. of res. b. 77b, 1973, pp 125-132.
!
!
!      latest modification: may 18, 1982
!
!
  implicit none

  integer ize, l, magx, n, nb, nbmx, ncalc, nend, nsig, nstart
  double precision alpha, b, em, empal, emp2al, en, enmten, ensig
  double precision enten, exparg, dgamma, half, halfx, one, p, plast, pold, psave
  double precision psavel, rtnsig, sum, tempa, tempb, tempc, test, tover, two, x
  double precision xlarge, zero
      dimension b(nb)
!
!  mathematical constants
!
      data one, two, zero, half /1.0d0,2.0d0,0.0d0,0.5d0/
!
!  machine dependent parameters
!
      data nsig, xlarge, exparg /17,1.0d4,88.0d0/
      data enten, ensig, rtnsig /1.0d38,1.0d17,1.0d-4/
      data enmten /1.2d-37/

      magx = ifix(sngl(x))
!
!  error return -- x,nb,or ize is out of range
!
      if ( &
        ( nb < 0 ) .or. &
        ( x < zero) .or. &
        ( alpha < zero) .or. &
        ( one <= alpha ) .or. &
        ( ize == 1 .and. exparg < x ) .or.  &
        ( ize == 2 .and. xlarge < x ) ) then
        ncalc = min ( nb, 0 ) - 1
        return
      end if
!
! use 2-term ascending series for small x
!
      ncalc = nb

      if (x.lt.rtnsig) go to 210
!
! initialize the forward sweep, the p-sequence of olver
!
      nbmx = nb - magx
      n = magx + 1
      en = dble(float(n+n)) + (alpha+alpha)
      plast = one
      p = en/x
!
! calculate general significance test
!
      test = ensig + ensig
      if ( 5 * nsig < 2 * magx ) then
        test = sqrt ( test * p )
      else
        test = test / 1.585d0**magx
      end if
!
! calculate p-sequence until n = nb-1.  check for possible overflow.
!
      if ( 3 <= nbmx ) then

        tover = enten/ensig
        nstart = magx + 2
        nend = nb - 1
        do n=nstart,nend
          en = en + two
          pold = plast
          plast = p
          p = en*plast/x + pold
          if (p.gt.tover) go to 40
        end do
        n = nend
        en = dble( n+n ) + (alpha+alpha)
!
! calculate special significance test for nbmx.gt.2.
!
        test = dmax1(test,dsqrt(plast*ensig)*dsqrt(p+p))

      end if
!
! calculate p-sequence until significance test passes
!
      do
        n = n + 1
        en = en + two
        pold = plast
        plast = p
        p = en * plast / x + pold
        if ( test <= p ) then
          exit
        end if
      end do

      go to 80
!
! to avoid overflow, divide p-sequence by tover.  calculate
! p-sequence until abs(p).gt.1.
!
   40 continue

      tover = enten
      p = p/tover
      plast = plast/tover
      psave = p
      psavel = plast
      nstart = n + 1

      do

        n = n + 1
        en = en + two
        pold = plast
        plast = p
        p = en*plast/x + pold

        if ( one < p ) then
          exit
        end if

      end do
!
! calculate backward test, and find ncalc, the highest n
! such that the test is passed.
!
      tempb = en/x
      test = pold*plast*(half-half/(tempb*tempb))/ensig
      p = plast*tover
      n = n - 1
      en = en - two
      nend = min(nb,n)
      do l=nstart,nend
        ncalc = l
        pold = psavel
        psavel = psave
        psave = en*psavel/x + pold
        if (psave*psavel.gt.test) go to 70
      end do
      ncalc = nend + 1
   70 ncalc = ncalc - 1
!
! initialize the backward recursion and the normalization sum
!
   80 continue

      n = n + 1
      en = en + two
      tempb = zero
      tempa = one/p
      em = dble ( n ) - one
      empal = em + alpha
      emp2al = (em-one) + (alpha+alpha)
      sum = tempa*empal*emp2al/em
      nend = n - nb

      if ( nend < 0 ) then
        go to 130
      end if

      if ( nend == 0 ) then
        go to 110
      end if
!
! recur backward via difference equation, calculating (but
! not storing) b(n), until n = nb.
!
      do l=1,nend
        n = n - 1
        en = en - two
        tempc = tempb
        tempb = tempa
        tempa = (en*tempb)/x + tempc
        em = em - one
        emp2al = emp2al - one
        if ( n == 1 ) then
          exit
        end if
        if (n == 2) emp2al = one
        empal = empal - one
        sum = (sum+tempa*empal)*emp2al/em
      end do
!
! store b(nb)
!
  110 continue

      b(n) = tempa
      if (nb.gt.1) go to 120
      sum = (sum+sum) + tempa
      go to 190
!
! calculate and store b(nb-1)
!
  120 continue

      n = n - 1
      en = en - two
      b(n) = (en*tempa)/x + tempb
      if (n == 1) go to 180
      em = em - one
      emp2al = emp2al - one
      if (n == 2) emp2al = one
      empal = empal - one
      sum = (sum+b(n)*empal)*emp2al/em
      go to 150
!
! n.lt.nb, so store b(n) and set higher orders to zero
!
  130 continue

      b(n) = tempa
      nend = -nend
      do l=1,nend
        b(n+l) = zero
      end do

  150 continue

      nend = n - 2
!
! calculate via difference equation and store b(n), until n = 2
!
      do l=1,nend
        n = n - 1
        en = en - two
        b(n) = (en*b(n+1))/x + b(n+2)
        em = em - one
        emp2al = emp2al - one
        if (n == 2) emp2al = one
        empal = empal - one
        sum = (sum+b(n)*empal)*emp2al/em
      end do
!
! calculate b(1)
!
      b(1) = two*empal*b(2)/x + b(3)

  180 continue

      sum = (sum+sum) + b(1)
!
! normalize.  divide all b(n) by sum.
!
  190 continue

      if (alpha.ne.zero) then
        sum = sum*dgamma(one+alpha)*(x*half)**(-alpha)
      end if

      if (ize == 1) then
        sum = sum*dexp(-x)
      end if

      tempa = enmten
      if (sum.gt.one) then
        tempa = tempa*sum
      end if

      do n = 1, nb
        if ( b(n) < tempa ) then
          b(n) = zero
        else
          b(n) = b(n) / sum
        end if
      end do

      return
!
! two-term ascending series for small x
!
  210 continue

      tempa = one
      empal = one + alpha
      halfx = zero
      if (x.gt.enmten) halfx = half*x
      if (alpha.ne.zero) tempa = halfx**alpha/dgamma(empal)
      if (ize == 2) tempa = tempa*dexp(-x)
      tempb = zero
      if ((x+one).gt.one) tempb = halfx*halfx
      b(1) = tempa + tempa*tempb/empal

      if ((x.ne.zero) .and. (b(1) == zero)) then
        ncalc = 0
      end if

      if ( nb == 1 ) then
        return
      end if

      if ( x <= zero ) then

        b(2:nb) = zero
!
! calculate higher order functions
!
      else

        tempc = halfx
        tover = ( enmten + enmten ) / x

        if ( tempb .ne. zero ) then
          tover = enmten / tempb
        end if

        do n = 2, nb

          tempa = tempa / empal
          empal = empal + one
          tempa = tempa * tempc
          if ( tempa .le. tover * empal ) then
            tempa = zero
          end if
          b(n) = tempa + tempa * tempb / empal

          if ( ( b(n) == zero ) .and. ( ncalc .gt. n ) ) then
            ncalc = n - 1
          end if

        end do

      end if

      return
      end