function clausen ( x )

c*********************************************************************72
c
cc clausen() evaluates the Clausen function Cl2(x).
c
c  Discussion:
c
c    Note that the first coefficient, a0 in Koelbig's paper,
c    is doubled here, to account for a different convention in
c    Chebyshev coefficients.
c
c    The program was producing NaN for x = 0 or multiples of 2 pi.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    30 March 2025
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Kurt Koelbig,
c    Chebyshev coefficients for the Clausen function Cl2(x),
c    Journal of Computational and Applied Mathematics,
c    Volume 64, Number 3, 1995, pages 295-297.
c
c  Input:
c
c    double precision X: the evaluation point.
c
c  Output:
c
c    double precision CLAUSEN: the value of the function.
c
      implicit none
c
c  Chebyshev expansion for -pi/2 < x < +pi/2.
c
      double precision :: c1(19) = (/
     &  0.05590566394715132269D+00,
     &  0.00000000000000000000D+00,
     &  0.00017630887438981157D+00,  
     &  0.00000000000000000000D+00,
     &  0.00000126627414611565D+00,    
     &  0.00000000000000000000D+00,
     &  0.00000001171718181344D+00,
     &  0.00000000000000000000D+00,
     &  0.00000000012300641288D+00,
     &  0.00000000000000000000D+00,
     &  0.00000000000139527290D+00,
     &  0.00000000000000000000D+00,
     &  0.00000000000001669078D+00,
     &  0.00000000000000000000D+00,
     &  0.00000000000000020761D+00,
     &  0.00000000000000000000D+00,
     &  0.00000000000000000266D+00,
     &  0.00000000000000000000D+00,
     &  0.00000000000000000003D+00 /)
c
c  Chebyshev expansion for pi/2 < x < 3pi/2.
c
      double precision :: c2(32) = (/
     &   0.00000000000000000000D+00, 
     &  -0.96070972149008358753D+00, 
     &   0.00000000000000000000D+00, 
     &   0.04393661151911392781D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00078014905905217505D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00002621984893260601D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00000109292497472610D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00000005122618343931D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00000000258863512670D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00000000013787545462D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00000000000763448721D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00000000000043556938D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00000000000002544696D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00000000000000151561D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00000000000000009172D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00000000000000000563D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00000000000000000035D+00, 
     &   0.00000000000000000000D+00, 
     &   0.00000000000000000002D+00 /)

      double precision clausen
      integer, parameter :: n1 = 19
      integer, parameter :: n2 = 30
      double precision r8_csevl
      double precision, parameter :: r8_pi = 3.141592653589793D+00
      double precision value
      double precision x
      double precision x2
      double precision x3
      double precision xa
      double precision xb
      double precision xc
c
c  The function is periodic.  Wrap X into [-pi/2, 3pi/2].
c
      xa = - 0.5 * r8_pi
      xb =   0.5 * r8_pi
      xc =   1.5 * r8_pi

      x2 = x
      do while ( x2 < xa )
        x2 = x2 + 2.0D+00 * r8_pi
      end do

      do while ( xc < x2 )
        x2 = x2 - 2.0D+00 * r8_pi
      end do
c
c  Choose the appropriate expansion.
c
      if ( abs ( x2 ) < epsilon ( x2 ) ) then
        value = 0.0D+00
      else if ( x2 < xb ) then
        x3 = 2.0D+00 * x2 / r8_pi
        value = x2 - x2 * log ( abs ( x2 ) )
     &    + 0.5D+00 * x2 ** 3 * r8_csevl ( x3, c1, n1 )
      else
        x3 = 2.0D+00 * x2 / r8_pi - 2.0D+00
        value = r8_csevl ( x3, c2, n2 )
      end if

      clausen = value

      return
      end
      function r8_csevl ( x, a, n )

c*********************************************************************72
c
cc r8_csevl() evaluates a Chebyshev series.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    06 September 2021
c
c  Author:
c
c    John Burkardt.
c
c  Reference:
c
c    Roger Broucke,
c    Algorithm 446:
c    Ten Subroutines for the Manipulation of Chebyshev Series,
c    Communications of the ACM,
c    Volume 16, Number 4, April 1973, pages 254-256.
c
c  Input:
c
c    double precision X, the evaluation point.
c
c    double precision A(N), the Chebyshev coefficients.
c
c    integer N, the number of Chebyshev coefficients.
c
c  Output:
c
c    double precision R8_CSEVL, the Chebyshev series evaluated at X.
c
      implicit none

      integer n

      double precision a(n)
      double precision b0
      double precision b1
      double precision b2
      integer i
      double precision r8_csevl
      double precision twox
      double precision x

      if ( n < 1 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'R8_CSEVL - Fatal error!'
        write ( *, '(a)' ) '  Number of terms <= 0.'
        stop 1
      end if

      if ( 1000 < n ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'R8_CSEVL - Fatal error!'
        write ( *, '(a)' ) '  1000 < Number of terms.'
        stop 1
      end if

      if ( x < -1.1D+00 .or. 1.1D+00 < x ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'R8_CSEVL - Fatal error!'
        write ( *, '(a)' ) '  X outside [-1,+1]'
        write ( *, '(a,g14.6)' ) '  X = ', x
        stop 1
      end if

      twox = 2.0D+00 * x
      b1 = 0.0D+00
      b0 = 0.0D+00

      do i = n, 1, -1
        b2 = b1
        b1 = b0
        b0 = twox * b1 - b2 + a(i)
      end do

      r8_csevl = 0.5D+00 * ( b0 - b2 )

      return
      end