function si ( x )

c*********************************************************************72
c
cc si() computes the sine integral 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    09 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 si: the value of Si(x).
c
      implicit none

      double precision bj(101)
      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 xabs
      double precision xcs
      double precision xf
      double precision xg
      double precision xg1
      double precision xg2
      double precision xr
      double precision xs
      double precision xsign
      double precision xss

      p2 = 1.570796326794897D+00
      el = 0.5772156649015329D+00
      eps = 1.0D-15
      x2 = x * x
      xabs = abs ( x )
      if ( x < 0.0 ) then
        xsign = -1.0D+00
      else
        xsign = +1.0D+00
      end if

      if ( xabs == 0.0D+00 ) then

        si = 0.0D+00

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

        xr = xabs
        si = xabs
        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
            si = xsign * si
            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 )
        si = xsign * ( xabs * xcs * xg1 
     &    + 2.0 * xss * xg2 - sin ( xabs ) )

      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
        si = xsign * ( p2 - xf * cos ( xabs ) / xabs 
     &    - xg * sin ( xabs ) / xabs )

      end if

      return
      end