function value = clausen ( x )

%*****************************************************************************80
%
%% CLAUSEN evaluates the Clausen function Cl2(x).
%
%  Discussion:
%
%    Note that the first coefficient, a0 in Koelbig's paper, 
%    is doubled here, to account for a different convention in
%    Chebyshev coefficients.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 December 2016
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Kurt Koelbig,
%    Chebyshev coefficients for the Clausen function Cl2(x),
%    Journal of Computational and Applied Mathematics,
%    Volume 64, Number 3, 1995, pages 295-297.
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of the function.
%

%
%  Chebyshev expansion for -pi/2 < x < +pi/2.
%
  n1 = 19;
  c1 = [ ...
    0.05590566394715132269; ...
    0.00000000000000000000; ...
    0.00017630887438981157; ...
    0.00000000000000000000; ...
    0.00000126627414611565; ...
    0.00000000000000000000; ...
    0.00000001171718181344; ...
    0.00000000000000000000; ...
    0.00000000012300641288; ...
    0.00000000000000000000; ...
    0.00000000000139527290; ...
    0.00000000000000000000; ...
    0.00000000000001669078; ...
    0.00000000000000000000; ...
    0.00000000000000020761; ...
    0.00000000000000000000; ...
    0.00000000000000000266; ...
    0.00000000000000000000; ...
    0.00000000000000000003 ];
%
%  Chebyshev expansion for pi/2 < x < 3pi/2.
%
  n2 = 32;
  c2 = [ ...
    0.00000000000000000000; ...
   -0.96070972149008358753; ...
    0.00000000000000000000; ...
    0.04393661151911392781; ...
    0.00000000000000000000; ...
    0.00078014905905217505; ...
    0.00000000000000000000; ...
    0.00002621984893260601; ...
    0.00000000000000000000; ...
    0.00000109292497472610; ...
    0.00000000000000000000; ...
    0.00000005122618343931; ...
    0.00000000000000000000; ...
    0.00000000258863512670; ...
    0.00000000000000000000; ...
    0.00000000013787545462; ...
    0.00000000000000000000; ...
    0.00000000000763448721; ...
    0.00000000000000000000; ...
    0.00000000000043556938; ...
    0.00000000000000000000; ...
    0.00000000000002544696; ...
    0.00000000000000000000; ...
    0.00000000000000151561; ...
    0.00000000000000000000; ...
    0.00000000000000009172; ...
    0.00000000000000000000; ...
    0.00000000000000000563; ...
    0.00000000000000000000; ...
    0.00000000000000000035; ...
    0.00000000000000000000; ...
    0.00000000000000000002 ];
%
%  The function is periodic.  Wrap X into [-pi/2, 3pi/2].
%
  xa = - 0.5 * pi;
  xb =   0.5 * pi;
  xc =   1.5 * pi;
  x2 = r8_wrap ( x, xa, xc );
%
%  Choose the appropriate expansion.
%
  if ( x2 < xb )
    x3 = 2.0 * x2 / pi;
    value = x2 - x2 * log ( abs ( x2 ) ) + 0.5 * x2 ^ 3 * r8_csevl ( x3, c1, n1 );
  else
    x3 = 2.0 * x2 / pi - 2.0;
    value = r8_csevl ( x3, c2, n2 );
  end

  return
end

