ik01b Subroutine

subroutine ik01b ( x, bi0, di0, bi1, di1, bk0, dk0, bk1, dk1 )

  !*****************************************************************************80
  !
  !! IK01B: Bessel functions I0(x), I1(x), K0(x), and K1(x) and derivatives.
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    17 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) BI0, DI0, BI1, DI1, BK0, DK0, BK1, DK1, the
  !    values of I0(x), I0'(x), I1(x), I1'(x), K0(x), K0'(x), K1(x), K1'(x).
  !
  implicit none

  real ( kind = 8 ) bi0
  real ( kind = 8 ) bi1
  real ( kind = 8 ) bk0
  real ( kind = 8 ) bk1
  real ( kind = 8 ) di0
  real ( kind = 8 ) di1
  real ( kind = 8 ) dk0
  real ( kind = 8 ) dk1
  real ( kind = 8 ) t
  real ( kind = 8 ) t2
  real ( kind = 8 ) x

  if ( x == 0.0D+00 ) then

     bi0 = 1.0D+00
     bi1 = 0.0D+00
     bk0 = 1.0D+300
     bk1 = 1.0D+300
     di0 = 0.0D+00
     di1 = 0.5D+00
     dk0 = -1.0D+300
     dk1 = -1.0D+300
     return

  else if ( x <= 3.75D+00 ) then

     t = x / 3.75D+00
     t2 = t * t

     bi0 = ((((( &
          0.0045813D+00   * t2 &
          + 0.0360768D+00 ) * t2 &
          + 0.2659732D+00 ) * t2 &
          + 1.2067492D+00 ) * t2 &
          + 3.0899424D+00 ) * t2 &
          + 3.5156229D+00 ) * t2 &
          + 1.0D+00

     bi1 = x * (((((( &
          0.00032411D+00   * t2 &
          + 0.00301532D+00 ) * t2 &
          + 0.02658733D+00 ) * t2 &
          + 0.15084934D+00 ) * t2 &
          + 0.51498869D+00 ) * t2 &
          + 0.87890594D+00 ) * t2 &
          + 0.5D+00 )

  else

     t = 3.75D+00 / x

     bi0 = (((((((( &
          0.00392377D+00   * t &
          - 0.01647633D+00 ) * t &
          + 0.02635537D+00 ) * t &
          - 0.02057706D+00 ) * t &
          + 0.916281D-02 ) * t &
          - 0.157565D-02 ) * t &
          + 0.225319D-02 ) * t &
          + 0.01328592D+00 ) * t &
          + 0.39894228D+00 ) * exp ( x ) / sqrt ( x )

     bi1 = (((((((( &
          - 0.420059D-02     * t &
          + 0.01787654D+00 ) * t &
          - 0.02895312D+00 ) * t &
          + 0.02282967D+00 ) * t &
          - 0.01031555D+00 ) * t &
          + 0.163801D-02 ) * t &
          - 0.00362018D+00 ) * t &
          - 0.03988024D+00 ) * t &
          + 0.39894228D+00 ) * exp ( x ) / sqrt ( x )

  end if

  if ( x <= 2.0D+00 ) then

     t = x / 2.0D+00
     t2 = t * t

     bk0 = ((((( &
          0.0000074D+00   * t2 &
          + 0.0001075D+00 ) * t2 &
          + 0.00262698D+00 ) * t2 &
          + 0.0348859D+00 ) * t2 &
          + 0.23069756D+00 ) * t2 &
          + 0.4227842D+00 ) * t2 &
          - 0.57721566D+00 - bi0 * log ( t )

     bk1 = (((((( &
          - 0.00004686D+00   * t2 &
          - 0.00110404D+00 ) * t2 &
          - 0.01919402D+00 ) * t2 &
          - 0.18156897D+00 ) * t2 &
          - 0.67278579D+00 ) * t2 &
          + 0.15443144D+00 ) * t2 &
          + 1.0D+00 ) / x + bi1 * log ( t )

  else

     t = 2.0D+00 / x
     t2 = t * t

     bk0 = (((((( &
          0.00053208D+00   * t &
          - 0.0025154D+00 )  * t &
          + 0.00587872D+00 ) * t &
          - 0.01062446D+00 ) * t &
          + 0.02189568D+00 ) * t &
          - 0.07832358D+00 ) * t &
          + 1.25331414D+00 ) * exp ( - x ) / sqrt ( x )

     bk1 = (((((( &
          - 0.00068245D+00   * t &
          + 0.00325614D+00 ) * t &
          - 0.00780353D+00 ) * t &
          + 0.01504268D+00 ) * t &
          - 0.0365562D+00  ) * t & 
          + 0.23498619D+00 ) * t &
          + 1.25331414D+00 ) * exp ( - x ) / sqrt ( x )

  end if

  di0 = bi1
  di1 = bi0 - bi1 / x
  dk0 = -bk1
  dk1 = -bk0 - bk1 / x

  return
end subroutine ik01b