/***********************************************************************
* This program computes the value of the integral of func(x) over the  *
* interval (a,b) by using the summed Clenshaw-Curtis formula           *
* -------------------------------------------------------------------- *
* SAMPLE RUN:                                                          *
* (Integrate function func(x) = cos(x) - x*sin(x) from x=0 to x=PI)    *
*                                                                      *
* Here parameters m=4 (number of subintervals in [0, PI]),             *
*                 n=6 (number of nodes for Clenshaw-Curtis quadrature) *
*                                                                      *
* Integral = -3.141593                                                 *
* Error code = 0                                                       *
*                                                                      *
* (Exact value is -PI).                                                *
* -------------------------------------------------------------------- *
* REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges        *
*        and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].           *
*                                                                      *
*                           Visual C++ Release By J-P Moreau, Paris.   *
*                                      (www.jpmoreau.fr)               *
***********************************************************************/
#include <stdio.h>
#include <math.h>

#define ONE   1.0
#define PI    4.0*atan(1.0)
#define REAL  double
#define ZERO  0.0

#define SIZE  10

//Function to integrate
REAL func(REAL x) {
  return (cos(x)-x*sin(x));
}

int ClenCurtStGew (
                   int n,
                   REAL* StStelle,
                   REAL* Gewicht
                  )
/***********************************************************************
* Computes the nodes and weights of a Clenshaw-Curtis quadrature       *
* formula of local error order  n+3 for the reference interval [-1,1]. *
*                                                                      *
* Parameters:                                                          *
*   int n                n+1 is the number of nodes and weights:       *
*                        n > 1, n odd                                  *
*   REAL   StStelle []   computed nodes                                *
*   REAL   Gewicht  []   computed weights                              *
*                                                                      *
* Return value :                                                       *
*   0:                   all ok                                        *
*   1:                   n too small or odd                            *
*                                                                      *
* Author:                Uli Eggermann,  9.24.1990                     *
***********************************************************************/
{
  int k, j, m;
  REAL p, f, g, h, i = ONE, d;
  if (n < 2 || n % 2 != 0) return 1;              /* n - Test         */
  m = n / 2;                                      /* m = n/2          */
  d = (REAL)(n * n - 1);                          /* d = n*n+1        */
  h = 2. / (n * d);                               /* h = 2/(n(n*n+1)) */
  f = 4. / n;                                     /* f = 4/n          */
  p = f * atan (1.);                              /* p = pi/n         */

  Gewicht  [0] =       Gewicht  [n] = 1 / d;       /*  left end point */
  StStelle [0] = -ONE; StStelle [n] = ONE;         /* right end point */

  for (k = 1; k <= m; k++, i = -i)                /* Interval (-1,0]  */
  {
    for (g = ZERO, j = 1; j < m; j++)
         g += cos(2 * j * k * p) / (4 * j * j - 1);
    Gewicht  [k] = h * (d + i) - f * g;
    StStelle [k] = - cos(k * p);
  }
  StStelle [m] = ZERO;                          /* center of interval */

  for (k = m+1; k < n; k++)                       /* Interval (0, 1)  */
  {
    Gewicht  [k] =   Gewicht  [n - k];
    StStelle [k] = - StStelle [n - k];
  }
  return 0;
} // ClenCurtStGew

int ClenCurt (
              int   m,
              REAL* t,
              int   n,
              REAL* Tk,
              REAL* Ak,
              REAL* Resultat
             )
/***********************************************************************
* Computes the value of the integral of func (x) over the interval     *
* (a,b) with the partition                                             *
*              t: a = t[0] < t[1] < .. < t[m] = b                      *
* by using the summed  Clenshaw-Curtis formula.                        *
* This program uses precomputed [0..n] weight vectors and Chebyshev    *
* nodes.                                                               *
*                                                                      *
* Parameter:                                                           *
*   int    m              number of subintervals                       *
*   REAL   t []           partition                                    *
*   int    n              n + 1 = number of nodes, n > 1, n even       *
*                                   (n + 2 = global error order)       *
*   REAL   Tk []          Chebyshev nodes                              *
*   REAL   Ak []          Weights                                      *
*                           Tk and Ak must be made available before    *
*                           calling this function by the procedure     *
*                           ClenCurtStGew for example                  *
*   REAL   *Resultat      Compute integral value                       *
*                                                                      *
* Return value :                                                       *
*   0:                    o.k.                                         *
*   1:                    improper number of nodes                     *
*   2:                    improper number of sub intervals             *
*                                                                      *
* Author:                 Uli Eggermann, 10.3.1991                     *
***********************************************************************/
{
  int    j, k;
  REAL   v, h, sum;

  if (n < 2 || n % 2 != 0) return 1;        /* n positive ? n even  ? */
  if (m < 1)               return 2;        /* partition              */

  for (*Resultat = ZERO, j = 0; j < m; j++) /* loop over intervals    */
  {
    v = 0.5 * (t [j+1] - t [j]);            /* half the interval size */
    h = 0.5 * (t [j+1] + t [j]);            /* Interval center        */

    for (sum = ZERO, k = 0; k <= n; k++)    /* Chebyshev loop         */
      sum += Ak [k] * func (v * Tk [k] + h);

    *Resultat += v * sum;
  }
  return 0;
}

void main() {

  int error, m=4, n=6;
  REAL t[SIZE], Tk[SIZE], Ak[SIZE];
  REAL Result;

  t[0]=ZERO; t[1]=PI/4.0; t[2]=PI/2.0; t[3]=3.0*PI/4.0; t[4]=PI;

  ClenCurtStGew(n,Tk,Ak);

  error = ClenCurt(m,t,n,Tk,Ak,&Result);

  printf("\n Integral = %f\n", Result);
  printf(" Error code = %d\n\n", error);

}

/* ------------------------ END clencurt.cpp ------------------------ */