itjya Subroutine

subroutine itjya ( x, tj, ty )

  !*****************************************************************************80
  !
  !! ITJYA computes integrals of Bessel functions J0(t) and Y0(t).
  !
  !  Discussion:
  !
  !    This procedure integrates Bessel functions J0(t) and Y0(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:
  !
  !    25 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 ) TJ, TY, the integrals of J0(t) and Y0(t) 
  !    from 0 to x.
  !
  implicit none

  real ( kind = 8 ) a(18)
  real ( kind = 8 ) a0
  real ( kind = 8 ) a1
  real ( kind = 8 ) af
  real ( kind = 8 ) bf
  real ( kind = 8 ) bg
  real ( kind = 8 ) el
  real ( kind = 8 ) eps
  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 ) tj
  real ( kind = 8 ) ty
  real ( kind = 8 ) ty1
  real ( kind = 8 ) ty2
  real ( kind = 8 ) x
  real ( kind = 8 ) x2
  real ( kind = 8 ) xp

  pi = 3.141592653589793D+00
  el = 0.5772156649015329D+00
  eps = 1.0D-12

  if ( x == 0.0D+00 ) then

     tj = 0.0D+00
     ty = 0.0D+00

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

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

     ty1 = ( el + log ( x / 2.0D+00 ) ) * tj
     rs = 0.0D+00
     ty2 = 1.0D+00
     r = 1.0D+00

     do k = 1, 60
        r = -0.25D+00 * r * ( 2 * k - 1.0D+00 ) / ( 2 * k + 1.0D+00 ) &
             / ( k * k ) * x2
        rs = rs + 1.0D+00 / k
        r2 = r * ( rs + 1.0D+00 / ( 2.0D+00 * k + 1.0D+00 ) )
        ty2 = ty2 + r2
        if ( abs ( r2 ) < abs ( ty2 ) * eps ) then
           exit
        end if
     end do

     ty = ( ty1 - x * ty2 ) * 2.0D+00 / pi

  else

     a0 = 1.0D+00
     a1 = 5.0D+00 / 8.0D+00
     a(1) = a1

     do k = 1, 16
        af = ( ( 1.5D+00 * ( k + 0.5D+00 ) * ( k + 5.0D+00 / 6.0D+00 ) &
             * a1 - 0.5D+00 * ( k + 0.5D+00 ) * ( k + 0.5D+00 )  &
             * ( k - 0.5D+00 ) * a0 ) ) / ( k + 1.0D+00 )
        a(k+1) = af
        a0 = a1
        a1 = af
     end do

     bf = 1.0D+00
     r = 1.0D+00
     do k = 1, 8
        r = -r / ( x * x )
        bf = bf + a(2*k) * r
     end do
     bg = a(1) / x
     r = 1.0D+00 / x
     do k = 1, 8
        r = -r / ( x * x )
        bg = bg + a(2*k+1) * r
     end do
     xp = x + 0.25D+00 * pi
     rc = sqrt ( 2.0D+00 / ( pi * x ) )
     tj = 1.0D+00 - rc * ( bf * cos ( xp ) + bg * sin ( xp ) )
     ty = rc * ( bg * cos ( xp ) - bf * sin ( xp ) )

  end if

  return
end subroutine itjya