ittjya Subroutine

subroutine ittjya ( x, ttj, tty )

  !*****************************************************************************80
  !
  !! ITTJYA integrates (1-J0(t))/t from 0 to x, and Y0(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 ) TTJ, TTY, the integrals of [1-J0(t)]/t 
  !    from 0 to x and of Y0(t)/t from x to oo.
  !
  implicit none

  real ( kind = 8 ) a0
  real ( kind = 8 ) b1
  real ( kind = 8 ) bj0
  real ( kind = 8 ) bj1
  real ( kind = 8 ) by0
  real ( kind = 8 ) by1
  real ( kind = 8 ) e0
  real ( kind = 8 ) el
  real ( kind = 8 ) g0
  real ( kind = 8 ) g1
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  real ( kind = 8 ) pi
  real ( kind = 8 ) px
  real ( kind = 8 ) qx
  real ( kind = 8 ) r
  real ( kind = 8 ) r0
  real ( kind = 8 ) r1
  real ( kind = 8 ) r2
  real ( kind = 8 ) rs
  real ( kind = 8 ) t
  real ( kind = 8 ) ttj
  real ( kind = 8 ) tty
  real ( kind = 8 ) vt
  real ( kind = 8 ) x
  real ( kind = 8 ) xk

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

  if ( x == 0.0D+00 ) then

     ttj = 0.0D+00
     tty = -1.0D+300

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

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

     ttj = ttj * 0.125D+00 * x * x
     e0 = 0.5D+00 * ( pi * pi / 6.0D+00 - el * el ) &
          - ( 0.5D+00 * log ( x / 2.0D+00 ) + el ) &
          * log ( x / 2.0D+00 )
     b1 = el + log ( x / 2.0D+00 ) - 1.5D+00
     rs = 1.0D+00
     r = -1.0D+00
     do k = 2, 100
        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 ) < abs ( b1 ) * 1.0D-12 ) then
           exit
        end if
     end do

     tty = 2.0D+00 / pi * ( e0 + 0.125D+00 * x * x * b1 )

  else

     a0 = sqrt ( 2.0D+00 / ( pi * x ) )

     do l = 0, 1

        vt = 4.0D+00 * l * l
        px = 1.0D+00
        r = 1.0D+00
        do k = 1, 14
           r = - 0.0078125D+00 * r &
                * ( vt - ( 4.0D+00 * k - 3.0D+00 ) ** 2 ) &
                / ( x * k ) * ( vt - ( 4.0D+00 * k - 1.0D+00 ) ** 2 ) &
                / ( ( 2.0D+00 * k - 1.0D+00 ) * x )
           px = px + r
           if ( abs ( r ) < abs ( px ) * 1.0D-12 ) then
              exit
           end if
        end do

        qx = 1.0D+00
        r = 1.0D+00
        do k = 1, 14
           r = -0.0078125D+00 * r &
                * ( vt - ( 4.0D+00 * k - 1.0D+00 ) ** 2 ) &
                / ( x * k ) * ( vt - ( 4.0D+00 * k + 1.0D+00 ) ** 2 ) &
                / ( 2.0D+00 * k + 1.0D+00 ) / x
           qx = qx + r
           if ( abs ( r ) < abs ( qx ) * 1.0D-12 ) then
              exit
           end if
        end do

        qx = 0.125D+00 * ( vt - 1.0D+00 ) / x * qx
        xk = x - ( 0.25D+00 + 0.5D+00 * l ) * pi
        bj1 = a0 * ( px * cos ( xk ) - qx * sin ( xk ) )
        by1 = a0 * ( px * sin ( xk ) + qx * cos ( xk ) )
        if ( l == 0 ) then
           bj0 = bj1
           by0 = by1
        end if

     end do

     t = 2.0D+00 / x
     g0 = 1.0D+00
     r0 = 1.0D+00
     do k = 1, 10
        r0 = - k * k * t * t *r0
        g0 = g0 + r0
     end do

     g1 = 1.0D+00
     r1 = 1.0D+00
     do k = 1, 10
        r1 = - k * ( k + 1.0D+00 ) * t * t * r1
        g1 = g1 + r1
     end do

     ttj = 2.0D+00 * g1 * bj0 / ( x * x ) - g0 * bj1 / x &
          + el + log ( x / 2.0D+00 )
     tty = 2.0D+00 * g1 * by0 / ( x * x ) - g0 * by1 / x

  end if

  return
end subroutine ittjya