itjyb Subroutine

subroutine itjyb ( x, tj, ty )

  !*****************************************************************************80
  !
  !! ITJYB 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 ) f0
  real ( kind = 8 ) g0
  real ( kind = 8 ) pi
  real ( kind = 8 ) t
  real ( kind = 8 ) tj
  real ( kind = 8 ) ty
  real ( kind = 8 ) x
  real ( kind = 8 ) x1
  real ( kind = 8 ) xt

  pi = 3.141592653589793D+00

  if ( x == 0.0D+00 ) then

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

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

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

     tj = ((((((( &
          - 0.133718D-03      * t &
          + 0.2362211D-02 )   * t &
          - 0.025791036D+00 ) * t &
          + 0.197492634D+00 ) * t &
          - 1.015860606D+00 ) * t &
          + 3.199997842D+00 ) * t &
          - 5.333333161D+00 ) * t &
          + 4.0D+00 ) * x1

     ty = (((((((( &
          0.13351D-04       * t &
          - 0.235002D-03 )    * t &
          + 0.3034322d-02 )   * t &
          - 0.029600855D+00 ) * t &
          + 0.203380298D+00 ) * t &
          - 0.904755062D+00 ) * t &
          + 2.287317974D+00 ) * t &
          - 2.567250468D+00 ) * t &
          + 1.076611469D+00 ) * x1

     ty = 2.0D+00 / pi * log ( x / 2.0D+00 ) * tj - ty

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

     xt = x - 0.25D+00 * pi
     t = 16.0D+00 / ( x * x )

     f0 = (((((( &
          0.1496119D-02     * t &
          - 0.739083D-02 )    * t &
          + 0.016236617D+00 ) * t &
          - 0.022007499D+00 ) * t &
          + 0.023644978D+00 ) * t &
          - 0.031280848D+00 ) * t &
          + 0.124611058D+00 ) * 4.0D+00 / x

     g0 = ((((( &
          0.1076103D-02     * t &
          - 0.5434851D-02 )   * t &
          + 0.01242264D+00 )  * t &
          - 0.018255209D+00 ) * t &
          + 0.023664841D+00 ) * t &
          - 0.049635633D+00 ) * t &
          + 0.79784879D+00

     tj = 1.0D+00 - ( f0 * cos ( xt ) - g0 * sin ( xt ) ) / sqrt ( x )

     ty = - ( f0 * sin ( xt ) + g0 * cos ( xt ) ) / sqrt ( x )

  else

     t = 64.0D+00 / ( x * x )
     xt = x-0.25D+00 * pi

     f0 = ((((((( &
          - 0.268482D-04     * t &
          + 0.1270039D-03 )  * t &
          - 0.2755037D-03 )  * t &
          + 0.3992825D-03 )  * t &
          - 0.5366169D-03 )  * t &
          + 0.10089872D-02 ) * t &
          - 0.40403539D-02 ) * t &
          + 0.0623347304D+00 ) * 8.0D+00 / x

     g0 = (((((( &
          - 0.226238D-04        * t &
          + 0.1107299D-03 )     * t &
          - 0.2543955D-03 )     * t &
          + 0.4100676D-03 )     * t &
          - 0.6740148D-03 )     * t &
          + 0.17870944D-02 )    * t &
          - 0.01256424405D+00 ) * t &
          + 0.79788456D+00

     tj = 1.0D+00  - ( f0 * cos ( xt ) - g0 * sin ( xt ) ) / sqrt ( x )

     ty = - ( f0 * sin ( xt ) + g0 * cos ( xt ) ) / sqrt ( x )

  end if

  return
end subroutine itjyb