subroutine cisi ( x, ci, si )

c*********************************************************************72
c
cc cisi() computes the cosine integral Ci(x) and 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    28 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
      function sincn_antideriv ( x )

c*********************************************************************72
c
cc sincn_antideriv(): antiderivative of the normalized sinc() function.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    28 March 2025
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    real x: the argument of the function.
c
c  Output:
c
c    real sincn_antideriv: the value of the antiderivative.
c
      implicit none

      double precision ci
      double precision, parameter :: r8_pi = 3.141592653589793D+00
      double precision si
      double precision sincn_antideriv
      double precision x

      call cisi ( r8_pi * x, ci, si )
      sincn_antideriv = si / r8_pi

      return
      end
      function sincn_deriv2 ( x )

c*********************************************************************72
c
cc sincn_deriv2(): second derivative of the normalized sinc() function.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    28 March 2025
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    real x: the argument of the function.
c
c  Output:
c
c    real sincn_deriv2: the value of the second derivative.
c
      implicit none

      double precision, parameter :: r8_pi = 3.141592653589793D+00
      double precision sincn_deriv2
      double precision x

      if ( abs ( x ) <= 1.0D-10 ) then
        sincn_deriv2 = - r8_pi**2 / 3.0D+00
      else
        sincn_deriv2 =
     &    ( ( 2.0 - r8_pi**2 * x**2 ) * sin ( r8_pi * x )
     &      - 2.0 * r8_pi    * x      * cos ( r8_pi * x ) ) 
     &    / ( r8_pi * x**3 )
      end if

      return
      end
      function sincn_deriv ( x )

c*********************************************************************72
c
cc sincn_deriv() evaluates the derivative of the normalized sinc() function.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    28 March 2025
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    real x: the argument of the function.
c
c  Output:
c
c    real sincn_deriv: the value of the derivative.
c
      implicit none

      double precision, parameter :: r8_pi = 3.141592653589793D+00
      double precision sincn_deriv
      double precision x

      if ( x == 0.0D+00 ) then
        sincn_deriv = 0.0
      else
        sincn_deriv = ( r8_pi * x * cos ( r8_pi * x )                   
     &    - sin ( r8_pi * x ) ) / ( r8_pi * x**2 )
      end if

      return
      end
      function sincn_fun ( x )

c*********************************************************************72
c
cc sincn_fun() evaluates the normalized sinc() function.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    28 March 2025
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    real x: the argument of the function.
c
c  Output:
c
c    real sincn_fun: the value of the function.
c
      implicit none

      double precision, parameter :: r8_pi = 3.141592653589793D+00
      double precision sincn_fun
      double precision x

      if ( x == 0.0 ) then
        sincn_fun = 1.0
      else
        sincn_fun = sin ( r8_pi * x ) / ( r8_pi * x )
      end if

      return
      end
      function sincu_antideriv ( x )

c*********************************************************************72
c
cc sincu_antideriv(): antiderivative of the unnormalized sinc() function.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    28 March 2025
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    real x: the argument of the function.
c
c  Output:
c
c    real sincu_antideriv: the value of the antiderivative.
c
      implicit none

      double precision ci
      double precision sincu_antideriv
      double precision si
      double precision x

      call cisi ( x, ci, si )
      sincu_antideriv = si

      return
      end
      function sincu_deriv2 ( x )

c*********************************************************************72
c
cc sincu_deriv2(): second derivative of the unnormalized sinc() function.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    28 March 2025
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    real x: the argument of the function.
c
c  Output:
c
c    real sincu_deriv2: the value of the second derivative.
c
      implicit none

      double precision sincu_deriv2
      double precision x

      if ( x == 0.0D+00 ) then
        sincu_deriv2 = - 1.0D+00 / 3.0D+00
      else
        sincu_deriv2 = ( ( 2.0D+00 - x**2 ) * sin ( x )
     &    - 2.0 * x * cos ( x ) ) / x**3
      end if

      return
      end
      function sincu_deriv ( x )

c*********************************************************************72
c
cc sincu_deriv() evaluates the derivative of the unnormalized sinc() function.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    28 March 2025
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    real x: the argument of the function.
c
c  Output:
c
c    real sincu_deriv: the value of the derivative.
c
      implicit none

      double precision sincu_deriv
      double precision x

      if ( x == 0.0D+00 ) then
        sincu_deriv = 0.0D+00
      else
        sincu_deriv = ( x * cos ( x ) - sin ( x ) ) / x**2
      end if

      return
      end
      function sincu_fun ( x )

c*********************************************************************72
c
cc sincu_fun() evaluates the unnormalized sinc() function.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    28 March 2025
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    real x: the argument of the function.
c
c  Output:
c
c    real sincu_fun: the value of the function.
c
      implicit none

      double precision sincu_fun
      double precision x

      if ( x == 0.0D+00 ) then
        sincu_fun = 1.0D+00
      else
        sincu_fun = sin ( x ) / x
      end if

      return
      end