Prev Next dbskes.f Headings

      subroutine dbskes(xnu,x,nin,bke)
c***begin prologue  dbskes
c***date written   770601   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  c10b3
c***keywords  bessel function,double precision,exponentially scaled,
c             fractional order,modified bessel function,sequence,
c             special function,third kind
c***author  fullerton, w., (lanl)
c***purpose  computes a d.p. sequence of exponentially scaled modified
c            bessel functions of the third kind of fractional order.
c dbskes(xnu,x,nin,bke) computes a double precision sequence
c of exponentially scaled modified bessel functions
c of the third kind of order xnu + i at x, where x .gt. 0,
c xnu lies in (-1,1), and i = 0, 1, ... , nin - 1, if nin is positive
c and i = 0, -1, ... , nin + 1, if nin is negative.  on return, the
c vector bke(.) contains the results at x for order starting at xnu.
c xnu, x, and bke are double precison.  nin is integer.
c***references  (none)
c***routines called  d1mach,d9knus,xerror
c***end prologue  dbskes
      double precision xnu, x, bke(1), bknu1, v, vincr, vend, alnbig,
     1  d1mach
      data alnbig / 0.d0 /
c***first executable statement  dbskes
      if (alnbig.eq.0.d0) alnbig = dlog (d1mach(2))
      v = dabs(xnu)
      n = iabs(nin)
      if ( call xerror ( 'dbskes  abs(xnu) must be lt 1', 29,
     1  2, 2)
      if (x.le.0.d0) call xerror ( 'dbskes  x is le 0', 17, 3, 2)
      if (n.eq.0) call xerror ( 'dbskes  n the number in the sequence is
     1 0', 41, 4, 2)
      call d9knus (v, x, bke(1), bknu1, iswtch)
      if (n.eq.1) return
      vincr = sign (1.0, float(nin))
      direct = vincr
      if ( direct = vincr*dsign(1.d0, xnu)
      if (iswtch.eq.1 .and. call xerror ( 'dbskes  x so sm
     1all bessel k-sub-xnu+1 overflows', 47, 5, 2)
      bke(2) = bknu1
      if ( call d9knus (dabs(xnu+vincr), x, bke(2), bknu1,
     1  iswtch)
      if (n.eq.2) return
      vend = dabs (xnu+dble(float(nin))) - 1.0d0
      if ((vend-.5d0)*dlog(vend)+0.27d0-vend*(dlog(x)-.694d0).gt.alnbig)
     1  call xerror (  'dbskes  x so small or abs(nu) so big that bessel
     2 k-sub-nu overflows', 67, 5, 2)
      v = xnu
      do 10 i=3,n
        v = v + vincr
        bke(i) = 2.0d0*v*bke(i-1)/x + bke(i-2)
 10   continue

Input File: f.omh/dbskes.f.omh