itika Subroutine

subroutine itika ( x, ti, tk )

  !*****************************************************************************80
  !
  !! ITIKA computes the integral of the modified Bessel functions I0(t) and K0(t).
  !
  !  Discussion:
  !
  !    This procedure integrates modified 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:
  !
  !    18 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 integrals of I0(t) and K0(t)
  !    from 0 to X.
  !
  implicit none

  real ( kind = 8 ), save, dimension ( 10 ) :: a = (/ &
       0.625D+00,           1.0078125D+00, &
       2.5927734375D+00,    9.1868591308594D+00, &
       4.1567974090576D+01, 2.2919635891914D+02, &
       1.491504060477D+03,  1.1192354495579D+04, &
       9.515939374212D+04,  9.0412425769041D+05 /)
  real ( kind = 8 ) b1
  real ( kind = 8 ) b2
  real ( kind = 8 ) e0
  real ( kind = 8 ) el
  integer ( kind = 4 ) k
  real ( kind = 8 ) pi
  real ( kind = 8 ) r
  real ( kind = 8 ) rc1
  real ( kind = 8 ) rc2
  real ( kind = 8 ) rs
  real ( kind = 8 ) ti
  real ( kind = 8 ) tk
  real ( kind = 8 ) tw
  real ( kind = 8 ) x
  real ( kind = 8 ) x2

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

  if ( x == 0.0D+00 ) then

     ti = 0.0D+00
     tk = 0.0D+00
     return

  else if ( x < 20.0D+00 ) then

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

     ti = ti * x

  else

     ti = 1.0D+00
     r = 1.0D+00
     do k = 1, 10
        r = r / x
        ti = ti + a(k) * r
     end do
     rc1 = 1.0D+00 / sqrt ( 2.0D+00 * pi * x )
     ti = rc1 * exp ( x ) * ti

  end if

  if ( x < 12.0D+00 ) then

     e0 = el + log ( x / 2.0D+00 )
     b1 = 1.0D+00 - e0
     b2 = 0.0D+00
     rs = 0.0D+00
     r = 1.0D+00
     do k = 1, 50
        r = 0.25D+00 * r * ( 2 * k - 1.0D+00 ) &
             / ( 2 * k + 1.0D+00 ) / ( k * k ) * x2
        b1 = b1 + r * ( 1.0D+00 / ( 2 * k + 1 ) - e0 )
        rs = rs + 1.0D+00 / k
        b2 = b2 + r * rs
        tk = b1 + b2
        if ( abs ( ( tk - tw ) / tk ) < 1.0D-12 ) then
           exit
        end if
        tw = tk
     end do

     tk = tk * x

  else

     tk = 1.0D+00
     r = 1.0D+00
     do k = 1, 10
        r = -r / x
        tk = tk + a(k) * r
     end do
     rc2 = sqrt ( pi / ( 2.0D+00 * x ) )
     tk = pi / 2.0D+00 - rc2 * tk * exp ( - x )

  end if

  return
end subroutine itika