ittikb Subroutine

subroutine ittikb ( x, tti, ttk )

  !*****************************************************************************80
  !
  !! ITTIKB integrates (I0(t)-1)/t from 0 to x, K0(t)/t from x to infinity.
  !
  !  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:
  !
  !    28 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 integral limit.
  !
  !    Output, real ( kind = 8 ) TTI, TTK, the integrals of
  !    [I0(t)-1]/t from 0 to x, and K0(t)/t from x to oo.
  !
  implicit none

  real ( kind = 8 ) e0
  real ( kind = 8 ) el
  real ( kind = 8 ) pi
  real ( kind = 8 ) t
  real ( kind = 8 ) t1
  real ( kind = 8 ) tti
  real ( kind = 8 ) ttk
  real ( kind = 8 ) x
  real ( kind = 8 ) x1

  pi = 3.141592653589793D+00
  el = 0.5772156649015329D+00

  if ( x == 0.0D+00 ) then

     tti = 0.0D+00

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

     x1 = x / 5.0D+00
     t = x1 * x1
     tti = ((((((( &
          0.1263D-03       * t &
          + 0.96442D-03 )    * t &
          + 0.968217D-02 )   * t &
          + 0.06615507D+00 ) * t &
          + 0.33116853D+00 ) * t &
          + 1.13027241D+00 ) * t &
          + 2.44140746D+00 ) * t &
          + 3.12499991D+00 ) * t

  else

     t = 5.0D+00 / x
     tti = ((((((((( &
          2.1945464D+00   * t &
          -  3.5195009D+00 ) * t &
          - 11.9094395D+00 ) * t &
          + 40.394734D+00  ) * t &
          - 48.0524115D+00 ) * t &
          + 28.1221478D+00 ) * t &
          -  8.6556013D+00 ) * t &
          +  1.4780044D+00 ) * t &
          -  0.0493843D+00 ) * t &
          +  0.1332055D+00 ) * t &
          +  0.3989314D+00
     tti = tti * exp ( x ) / ( sqrt ( x ) * x )

  end if

  if ( x == 0.0D+00 ) then

     ttk = 1.0D+300

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

     t1 = x / 2.0D+00
     t = t1 * t1
     ttk = ((((( &
          0.77D-06         * t &
          + 0.1544D-04 )     * t &
          + 0.48077D-03 )    * t &
          + 0.925821D-02 )   * t &
          + 0.10937537D+00 ) * t &
          + 0.74999993D+00 ) * t
     e0 = el + log ( x / 2.0D+00 )
     ttk = pi * pi / 24.0D+00 + e0 * ( 0.5D+00 * e0 + tti ) - ttk

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

     t = 2.0D+00 / x
     ttk = ((( &
          0.06084D+00    * t &
          - 0.280367D+00 ) * t &
          + 0.590944D+00 ) * t &
          - 0.850013D+00 ) * t &
          + 1.234684D+00
     ttk = ttk * exp ( - x ) / ( sqrt ( x ) * x )

  else

     t = 4.0D+00 / x
     ttk = ((((( &
          0.02724D+00     * t &
          - 0.1110396D+00 ) * t &
          + 0.2060126D+00 ) * t &
          - 0.2621446D+00 ) * t &
          + 0.3219184D+00 ) * t &
          - 0.5091339D+00 ) * t &
          + 1.2533141D+00
     ttk = ttk * exp ( - x ) / ( sqrt ( x ) * x )

  end if

  return
end subroutine ittikb