subroutine cisi ( x, ci, si )

c*********************************************************************72
c
cc cisi() computes cosine Ci(x) and sine integrals Si(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    26 March 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 of Ci(x) and Si(x).
c
c  Output:
c
c    double precision ci, si: the values of Ci(x) and Si(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 si
      double precision x
      double precision x2
      double precision xa
      double precision xa0
      double precision xa1
      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

      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