# include <float.h>
# include <math.h>
# include <stdio.h>
# include <stdlib.h>

# include "clausen.h"

/******************************************************************************/

double clausen ( double x )

/******************************************************************************/
/*
  Purpose:

    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 MIT license.

  Modified:

    12 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.

  Input:

    double X, the evaluation point.

  Output:

    double CLAUSEN, the value of the function.
*/
{
/*
  Chebyshev expansion for -pi/2 < x < +pi/2.
*/
  static double c1[19] = {
    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.
*/
  static double c2[32] = {
    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 };
  const double r8_pi = 3.141592653589793;
  int n1 = 19;
  int n2 = 30;
  double value;
  double x2;
  double x3;
  double xa;
  double xb;
  double xc;
/*
  The function is periodic.  Wrap X into [-pi/2, 3pi/2].
*/
  xa = - 0.5 * r8_pi;
  xb =   0.5 * r8_pi;
  xc =   1.5 * r8_pi;

  x2 = x;
  while ( x2 < xa )
  {
    x2 = x2 + 2.0 * r8_pi;
  }

  while ( xc < x2 )
  {
    x2 = x2 - 2.0 * r8_pi;
  }
/*
  Choose the appropriate expansion.
*/
  if ( fabs ( x2 ) < DBL_EPSILON )
  {
    value = 0.0;
  }
  else if ( x2 < xb )
  {
    x3 = 2.0 * x2 / r8_pi;
    value = x2 - x2 * log ( fabs ( x2 ) ) 
      + 0.5 * pow ( x2, 3 ) * r8_csevl ( x3, c1, n1 );
  }
  else
  {
    x3 = 2.0 * x2 / r8_pi - 2.0;
    value = r8_csevl ( x3, c2, n2 );
  }

  return value;
}
/******************************************************************************/

double r8_csevl ( double x, double a[], int n )

/******************************************************************************/
/*
  Purpose:

    r8_csevl() evaluates a Chebyshev series.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 January 2012

  Author:

    This version by John Burkardt.

  Reference:

    Roger Broucke,
    Algorithm 446:
    Ten Subroutines for the Manipulation of Chebyshev Series,
    Communications of the ACM,
    Volume 16, Number 4, April 1973, pages 254-256.

  Input:

    double X, the evaluation point.

    double CS[N], the Chebyshev coefficients.

    int N, the number of Chebyshev coefficients.

  Output:

    double R8_CSEVL, the Chebyshev series evaluated at X.
*/
{
  double b0;
  double b1;
  double b2;
  int i;
  double twox;
  double value;

  if ( n < 1 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_CSEVL - Fatal error!\n" );
    fprintf ( stderr, "  Number of terms <= 0.\n" );
    exit ( 1 );
  }

  if ( 1000 < n )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_CSEVL - Fatal error!\n" );
    fprintf ( stderr, "  Number of terms greater than 1000.\n" );
    exit ( 1 );
 }

  if ( x < -1.1 || 1.1 < x )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_CSEVL - Fatal error!\n" );
    fprintf ( stderr, "  X outside (-1,+1).\n" );
    exit ( 1 );
  }

  twox = 2.0 * x;
  b1 = 0.0;
  b0 = 0.0;

  for ( i = n - 1; 0 <= i; i-- )
  {
    b2 = b1;
    b1 = b0;
    b0 = twox * b1 - b2 + a[i];
  }

  value = 0.5 * ( b0 - b2 );

  return value;
}