cisia Subroutine

subroutine cisia ( x, ci, si )

  !*****************************************************************************80
  !
  !! CISIA computes cosine Ci(x) and sine integrals Si(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:
  !
  !    03 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 argument of Ci(x) and Si(x).
  !
  !    Output, real ( kind = 8 ) CI, SI, the values of Ci(x) and Si(x).
  !
  implicit none

  real ( kind = 8 ) bj(101)
  real ( kind = 8 ) ci
  real ( kind = 8 ) el
  real ( kind = 8 ) eps
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  real ( kind = 8 ) p2
  real ( kind = 8 ) si
  real ( kind = 8 ) x
  real ( kind = 8 ) x2
  real ( kind = 8 ) xa
  real ( kind = 8 ) xa0
  real ( kind = 8 ) xa1
  real ( kind = 8 ) xcs
  real ( kind = 8 ) xf
  real ( kind = 8 ) xg
  real ( kind = 8 ) xg1
  real ( kind = 8 ) xg2
  real ( kind = 8 ) xr
  real ( kind = 8 ) xs
  real ( kind = 8 ) xss

  p2 = 1.570796326794897D+00
  el = 0.5772156649015329D+00
  eps = 1.0D-15
  x2 = x * x

  if ( x == 0.0D+00 ) then

     ci = -1.0D+300
     si = 0.0D+00

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

     xr = -0.25D+00 * x2
     ci = el + log ( x ) + xr
     do k = 2, 40
        xr = -0.5D+00 * xr * ( k - 1 ) / ( k * k * ( 2 * k - 1 ) ) * x2
        ci = ci + xr
        if ( abs ( xr ) < abs ( ci ) * eps ) then
           exit
        end if
     end do

     xr = x
     si = x
     do k = 1, 40
        xr = -0.5D+00 * xr * ( 2 * k - 1 ) / k / ( 4 * k * k + 4 * k + 1 ) * x2
        si = si + xr
        if ( abs ( xr ) < abs ( si ) * eps ) then
           return
        end if
     end do

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

     m = int ( 47.2D+00 + 0.82D+00 * x )
     xa1 = 0.0D+00
     xa0 = 1.0D-100
     do k = m, 1, -1
        xa = 4.0D+00 * k * xa0 / x - xa1
        bj(k) = xa
        xa1 = xa0
        xa0 = xa
     end do
     xs = bj(1)
     do k = 3, m, 2
        xs = xs + 2.0D+00 * bj(k)
     end do
     bj(1) = bj(1) / xs
     do k = 2, m
        bj(k) = bj(k) / xs
     end do
     xr = 1.0D+00
     xg1 = bj(1)
     do k = 2, m
        xr = 0.25D+00 * xr * ( 2.0D+00 * k - 3.0D+00 ) **2 &
             / ( ( k - 1.0D+00 ) * ( 2.0D+00 * k - 1.0D+00 ) ** 2 ) * x
        xg1 = xg1 + bj(k) * xr
     end do

     xr = 1.0D+00
     xg2 = bj(1)
     do k = 2, m
        xr = 0.25D+00 * xr * ( 2.0D+00 * k - 5.0D+00 )**2 &
             / ( ( k-1.0D+00 ) * ( 2.0D+00 * k - 3.0D+00 ) ** 2 ) * x
        xg2 = xg2 + bj(k) * xr
     end do

     xcs = cos ( x / 2.0D+00 )
     xss = sin ( x / 2.0D+00 )
     ci = el + log ( x ) - x * xss * xg1 + 2.0 * xcs * xg2 - 2.0 * xcs * xcs
     si = x * xcs * xg1 + 2.0 * xss * xg2 - sin ( x )

  else

     xr = 1.0D+00
     xf = 1.0D+00
     do k = 1, 9
        xr = -2.0D+00 * xr * k * ( 2 * k - 1 ) / x2
        xf = xf + xr
     end do
     xr = 1.0D+00 / x
     xg = xr
     do k = 1, 8
        xr = -2.0D+00 * xr * ( 2 * k + 1 ) * k / x2
        xg = xg + xr
     end do
     ci = xf * sin ( x ) / x - xg * cos ( x ) / x
     si = p2 - xf * cos ( x ) / x - xg * sin ( x ) / x

  end if

  return
end subroutine cisia