# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <time.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.

  Parameters:

    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 = r8_wrap ( x, xa, xc );
/*
  Choose the appropriate expansion.
*/
  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;
}
/******************************************************************************/

void clausen_values ( int *n_data, double *x, double *fx )

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

    CLAUSEN_VALUES returns some values of the Clausen's integral.

  Discussion:

    The function is defined by:

      CLAUSEN(x) = Integral ( 0 <= t <= x ) -ln ( 2 * sin ( t / 2 ) ) dt

    The data was reported by McLeod.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    25 August 2004

  Author:

    John Burkardt

  Reference:

    Milton Abramowitz, Irene Stegun,
    Handbook of Mathematical Functions,
    National Bureau of Standards, 1964,
    ISBN: 0-486-61272-4,
    LC: QA47.A34.

    Allan McLeod,
    Algorithm 757:
    MISCFUN: A software package to compute uncommon special functions,
    ACM Transactions on Mathematical Software,
    Volume 22, Number 3, September 1996, pages 288-301.

  Parameters:

    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
    first call.  On each call, the routine increments N_DATA by 1, and
    returns the corresponding data; when there is no more data, the
    output value of N_DATA will be 0 again.

    Output, double *X, the argument of the function.

    Output, double *FX, the value of the function.
*/
{
# define N_MAX 20

  static double fx_vec[N_MAX] = {
      0.14137352886760576684E-01,
      0.13955467081981281934E+00,
     -0.38495732156574238507E+00,
      0.84831187770367927099E+00,
      0.10139591323607685043E+01,
     -0.93921859275409211003E+00,
      0.72714605086327924743E+00,
      0.43359820323553277936E+00,
     -0.98026209391301421161E-01,
     -0.56814394442986978080E+00,
     -0.70969701784448921625E+00,
      0.99282013254695671871E+00,
     -0.98127747477447367875E+00,
     -0.64078266570172320959E+00,
      0.86027963733231192456E+00,
      0.39071647608680211043E+00,
      0.47574793926539191502E+00,
      0.10105014481412878253E+01,
      0.96332089044363075154E+00,
     -0.61782699481929311757E+00 };

  static double x_vec[N_MAX] = {
       0.0019531250E+00,
       0.0312500000E+00,
      -0.1250000000E+00,
       0.5000000000E+00,
       1.0000000000E+00,
      -1.5000000000E+00,
       2.0000000000E+00,
       2.5000000000E+00,
      -3.0000000000E+00,
       4.0000000000E+00,
       4.2500000000E+00,
      -5.0000000000E+00,
       5.5000000000E+00,
       6.0000000000E+00,
       8.0000000000E+00,
     -10.0000000000E+00,
      15.0000000000E+00,
      20.0000000000E+00,
     -30.0000000000E+00,
      50.0000000000E+00 };

  if ( *n_data < 0 )
  {
    *n_data = 0;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *fx = 0.0;
  }
  else
  {
    *x  = x_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }

  return;
# undef N_MAX
}
/******************************************************************************/

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:

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

  Parameters:

    Input, double X, the evaluation point.

    Input, double CS[N], the Chebyshev coefficients.

    Input, 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;
}
/******************************************************************************/

double r8_wrap ( double r, double rlo, double rhi )

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

    R8_WRAP forces an R8 to lie between given limits by wrapping.

  Discussion:

    An R8 is a double value.

  Example:

    RLO = 4.0, RHI = 8.0

     R  Value

    -2     8
    -1     4
     0     5
     1     6
     2     7
     3     8
     4     4
     5     5
     6     6
     7     7
     8     8
     9     4
    10     5
    11     6
    12     7
    13     8
    14     4

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    04 July 2011

  Author:

    John Burkardt

  Parameters:

    Input, double R, a value.

    Input, double RLO, RHI, the desired bounds.

    Output, double R8_WRAP, a "wrapped" version of the value.
*/
{
  int n;
  double rhi2;
  double rlo2;
  double rwide;
  double value;
/*
  Guarantee RLO2 < RHI2.
*/
  if ( rlo <= rhi )
  {
    rlo2 = rlo;
    rhi2 = rhi;
  }
  else
  {
    rlo2 = rhi;
    rhi2 = rlo;
  }
/*
  Find the width.
*/
  rwide = rhi2 - rlo2;
/*
  Add enough copies of (RHI2-RLO2) to R so that the
  result ends up in the interval RLO2 - RHI2.
*/
  if ( rwide == 0.0 )
  {
    value = rlo;
  }
  else if ( r < rlo2 )
  {
    n = ( int ) ( ( rlo2 - r ) / rwide ) + 1;
    value = r + n * rwide;
    if ( value == rhi )
    {
      value = rlo;
    }
  }
  else
  {
    n = ( int ) ( ( r - rlo2 ) / rwide );
    value = r - n * rwide;
    if ( value == rlo )
    {
      value = rhi;
    }
  }
  return value;
}
/******************************************************************************/

void timestamp ( )

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

    TIMESTAMP prints the current YMDHMS date as a time stamp.

  Example:

    17 June 2014 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 June 2014

  Author:

    John Burkardt

  Parameters:

    None
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  printf ( "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}