ittjyb Subroutine

subroutine ittjyb ( x, ttj, tty )

  !*****************************************************************************80
  !
  !! ITTJYB 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:
  !
  !    01 August 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 ) e0
  real ( kind = 8 ) el
  real ( kind = 8 ) f0
  real ( kind = 8 ) g0
  real ( kind = 8 ) pi
  real ( kind = 8 ) t
  real ( kind = 8 ) t1
  real ( kind = 8 ) ttj
  real ( kind = 8 ) tty
  real ( kind = 8 ) x
  real ( kind = 8 ) x1
  real ( kind = 8 ) xt

  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 <= 4.0D+00 ) then

     x1 = x / 4.0D+00
     t = x1 * x1

     ttj = (((((( &
          0.35817D-04 * t &
          - 0.639765D-03 ) * t &
          + 0.7092535D-02 ) * t &
          - 0.055544803D+00 ) * t &
          + 0.296292677D+00 ) * t &
          - 0.999999326D+00 ) * t &
          + 1.999999936D+00 ) * t

     tty = ((((((( &
          - 0.3546D-05        * t &
          + 0.76217D-04 )     * t &
          - 0.1059499D-02 )   * t &
          + 0.010787555D+00 ) * t &
          - 0.07810271D+00 )  * t &
          + 0.377255736D+00 ) * t &
          - 1.114084491D+00 ) * t &
          + 1.909859297D+00 ) * t

     e0 = el + log ( x / 2.0D+00 )
     tty = pi / 6.0D+00 + e0 / pi * ( 2.0D+00 * ttj - e0 ) - tty

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

     xt = x + 0.25D+00 * pi
     t1 = 4.0D+00 / x
     t = t1 * t1

     f0 = ((((( &
          0.0145369D+00 * t &
          - 0.0666297D+00 ) * t &
          + 0.1341551D+00 ) * t &
          - 0.1647797D+00 ) * t &
          + 0.1608874D+00 ) * t &
          - 0.2021547D+00 ) * t &
          + 0.7977506D+00

     g0 = (((((( &
          0.0160672D+00   * t &
          - 0.0759339D+00 ) * t &
          + 0.1576116D+00 ) * t &
          - 0.1960154D+00 ) * t &
          + 0.1797457D+00 ) * t &
          - 0.1702778D+00 ) * t &
          + 0.3235819D+00 ) * t1

     ttj = ( f0 * cos ( xt ) + g0 * sin ( xt ) ) / ( sqrt ( x ) * x )
     ttj = ttj + el + log ( x / 2.0D+00 )
     tty = ( f0 * sin ( xt ) - g0 * cos ( xt ) ) / ( sqrt ( x ) * x )

  else

     t = 8.0D+00 / x
     xt = x + 0.25D+00 * pi

     f0 = ((((( &
          0.18118D-02    * t &
          - 0.91909D-02 )  * t &
          + 0.017033D+00 ) * t &
          - 0.9394D-03 )   * t &
          - 0.051445D+00 ) * t &
          - 0.11D-05 )     * t &
          + 0.7978846D+00

     g0 = ((((( &
          - 0.23731D-02     * t &
          + 0.59842D-02 )   * t &
          + 0.24437D-02 )   * t &
          - 0.0233178D+00 ) * t &
          + 0.595D-04 )     * t &
          + 0.1620695D+00 ) * t

     ttj = ( f0 * cos ( xt ) + g0 * sin ( xt ) )  &
          / ( sqrt ( x ) * x ) + el + log ( x / 2.0D+00 )
     tty = ( f0 * sin ( xt ) - g0 * cos ( xt ) )  &
          / ( sqrt ( x ) * x )

  end if

  return
end subroutine ittjyb