itikb Subroutine

subroutine itikb ( x, ti, tk )

  !*****************************************************************************80
  !
  !! ITIKB computes the integral of the Bessel functions I0(t) and K0(t).
  !
  !  Discussion:
  !
  !    This procedure integrates Bessel functions I0(t) and K0(t)
  !    with respect to t from 0 to x.
  !
  !  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:
  !
  !    24 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 upper limit of the integral.
  !
  !    Output, real ( kind = 8 ) TI, TK, the integral of I0(t) and K0(t)
  !    from 0 to X.
  !
  implicit none

  real ( kind = 8 ) pi
  real ( kind = 8 ) t
  real ( kind = 8 ) t1
  real ( kind = 8 ) ti
  real ( kind = 8 ) tk
  real ( kind = 8 ) x

  pi = 3.141592653589793D+00

  if ( x == 0.0D+00 ) then

     ti = 0.0D+00

  else if ( x < 5.0D+00 ) then

     t1 = x / 5.0D+00
     t = t1 * t1
     ti = (((((((( &
          0.59434D-03 * t &
          + 0.4500642D-02 ) * t &
          + 0.044686921D+00 ) * t &
          + 0.300704878D+00 ) * t &
          + 1.471860153D+00 ) * t &
          + 4.844024624D+00 ) * t &
          + 9.765629849D+00 ) * t &
          +10.416666367D+00 ) * t &
          + 5.0D+00 ) * t1

  else if ( 5.0D+00 <= x .and. x <= 8.0D+00 ) then

     t = 5.0D+00 / x
     ti = ((( &
          - 0.015166D+00 * t &
          - 0.0202292D+00 ) * t &
          + 0.1294122D+00 ) * t &
          - 0.0302912D+00 ) * t &
          + 0.4161224D+00
     ti = ti * exp ( x ) / sqrt ( x )

  else

     t = 8.0D+00 / x
     ti = ((((( &
          - 0.0073995D+00 * t &
          + 0.017744D+00 ) * t &
          - 0.0114858D+00 ) * t &
          + 0.55956D-02 ) * t &
          + 0.59191D-02 ) * t &
          + 0.0311734D+00 ) * t &
          + 0.3989423D+00
     ti = ti * exp ( x ) / sqrt ( x )

  end if

  if ( x == 0.0D+00 ) then

     tk = 0.0D+00

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

     t1 = x / 2.0D+00
     t = t1 * t1
     tk = (((((( &
          0.116D-05        * t &
          + 0.2069D-04 )     * t &
          + 0.62664D-03 )    * t &
          + 0.01110118D+00 ) * t &
          + 0.11227902D+00 ) * t &
          + 0.50407836D+00 ) * t &
          + 0.84556868D+00 ) * t1
     tk = tk - log ( x / 2.0D+00 ) * ti

  else if ( 2.0D+00 < x .and. x <= 4.0D+00 ) then

     t = 2.0D+00 / x
     tk = ((( &
          0.0160395D+00   * t &
          - 0.0781715D+00 ) * t &
          + 0.185984D+00 )  * t &
          - 0.3584641D+00 ) * t &
          + 1.2494934D+00
     tk = pi / 2.0D+00 - tk * exp ( - x ) / sqrt ( x )

  else if ( 4.0D+00 < x .and. x <= 7.0D+00 ) then

     t = 4.0D+00 / x
     tk = ((((( &
          0.37128D-02 * t &
          - 0.0158449D+00 ) * t &
          + 0.0320504D+00 ) * t &
          - 0.0481455D+00 ) * t &
          + 0.0787284D+00 ) * t &
          - 0.1958273D+00 ) * t &
          + 1.2533141D+00
     tk = pi / 2.0D+00 - tk * exp ( - x ) / sqrt ( x )

  else

     t = 7.0D+00 / x
     tk = ((((( &
          0.33934D-03      * t &
          - 0.163271D-02 )   * t &
          + 0.417454D-02 )   * t &
          - 0.933944D-02 )   * t &
          + 0.02576646D+00 ) * t &
          - 0.11190289D+00 ) * t &
          + 1.25331414D+00
     tk = pi / 2.0D+00 - tk * exp ( - x ) / sqrt ( x )

  end if

  return
end subroutine itikb