ittika Subroutine

subroutine ittika ( x, tti, ttk )

  !*****************************************************************************80
  !
  !! ITTIKA 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:
  !
  !    23 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 of K0(t)/t from x to oo.
  !
  implicit none

  real ( kind = 8 ) b1
  real ( kind = 8 ), save, dimension ( 8 ) :: c = (/ &
       1.625D+00, 4.1328125D+00, &
       1.45380859375D+01, 6.553353881835D+01, &
       3.6066157150269D+02, 2.3448727161884D+03, &
       1.7588273098916D+04, 1.4950639538279D+05 /)
  real ( kind = 8 ) e0
  real ( kind = 8 ) el
  integer ( kind = 4 ) k
  real ( kind = 8 ) pi
  real ( kind = 8 ) r
  real ( kind = 8 ) r2
  real ( kind = 8 ) rc
  real ( kind = 8 ) rs
  real ( kind = 8 ) tti
  real ( kind = 8 ) ttk
  real ( kind = 8 ) x

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

  if ( x == 0.0D+00 ) then
     tti = 0.0D+00
     ttk = 1.0D+300
     return
  end if

  if ( x < 40.0D+00 ) then
     tti = 1.0D+00
     r = 1.0D+00
     do k = 2, 50
        r = 0.25D+00 * r * ( k - 1.0D+00 ) / ( k * k * k ) * x * x
        tti = tti + r
        if ( abs ( r / tti ) < 1.0D-12 ) then
           exit
        end if
     end do

     tti = tti * 0.125D+00 * x * x

  else

     tti = 1.0D+00
     r = 1.0D+00
     do k = 1, 8
        r = r / x
        tti = tti + c(k) * r
     end do
     rc = x * sqrt ( 2.0D+00 * pi * x )
     tti = tti * exp ( x ) / rc

  end if

  if ( x <= 12.0D+00 ) then

     e0 = ( 0.5D+00 * log ( x / 2.0D+00 ) + el ) &
          * log ( x / 2.0D+00 ) + pi * pi / 24.0D+00 + 0.5D+00 * el * el
     b1 = 1.5D+00 - ( el + log ( x / 2.0D+00 ) )
     rs = 1.0D+00
     r = 1.0D+00
     do k = 2, 50
        r = 0.25D+00 * r * ( k - 1.0D+00 ) / ( k * k * k ) * x * x
        rs = rs + 1.0D+00 / k
        r2 = r * ( rs + 1.0D+00 / ( 2.0D+00 * k ) &
             - ( el + log ( x / 2.0D+00 ) ) )
        b1 = b1 + r2
        if ( abs ( r2 / b1 ) < 1.0D-12 ) then
           exit
        end if
     end do

     ttk = e0 - 0.125D+00 * x * x * b1

  else

     ttk = 1.0D+00
     r = 1.0D+00
     do k = 1, 8
        r = - r / x
        ttk = ttk + c(k) * r
     end do
     rc = x * sqrt ( 2.0D+00 / pi * x )
     ttk = ttk * exp ( - x ) / rc

  end if

  return
end subroutine ittika