function ci ( x )

c*********************************************************************72
c
cc ci() computes the cosine Ci(x).
c
c  Licensing:
c
c    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However,
c    they give permission to incorporate this routine into a user program
c    provided that the copyright is acknowledged.
c
c  Modified:
c
c    11 August 2025
c
c  Author:
c
c    Original FORTRAN77 version by Shanjie Zhang, Jianming Jin.
c    This version by John Burkardt.
c
c  Reference:
c
c    Shanjie Zhang, Jianming Jin,
c    Computation of Special Functions,
c    Wiley, 1996,
c    ISBN: 0-471-11963-6,
c    LC: QA351.C45.
c
c  Input:
c
c    double precision x: the argument.
c
c  Output:
c
c    double precision ci: the value of Ci(x).
c
      implicit none

      double precision bj(101)
      double precision ci
      double precision el
      double precision eps
      integer k
      integer m
      double precision p2
      double precision x
      double precision x2
      double precision xa
      double precision xa0
      double precision xa1
      double precision xabs
      double precision xcs
      double precision xf
      double precision xg
      double precision xg1
      double precision xg2
      double precision xr
      double precision xs
      double precision xss

      p2 = 1.570796326794897D+00
      el = 0.5772156649015329D+00
      eps = 1.0D-15
      x2 = x * x
      xabs = abs ( x )

      if ( xabs == 0.0D+00 ) then

        ci = -1.0D+300

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

        xr = -0.25D+00 * x2
        ci = el + log ( xabs ) + 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
            return
          end if
        end do

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

        m = int ( 47.2D+00 + 0.82D+00 * xabs )
        xa1 = 0.0D+00
        xa0 = 1.0D-100
        do k = m, 1, -1
          xa = 4.0D+00 * k * xa0 / xabs - 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 ) * xabs
          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 ) * xabs
          xg2 = xg2 + bj(k) * xr
        end do

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

      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 / xabs
        xg = xr
        do k = 1, 8
          xr = -2.0D+00 * xr * ( 2 * k + 1 ) * k / x2
          xg = xg + xr
        end do
        ci = xf * sin ( xabs ) / xabs - xg * cos ( xabs ) / xabs
 
      end if

      return
      end