DECLARE FUNCTION func# (x#)
DECLARE FUNCTION IClenCurtStGew% (n%, StStelle#(), Gewicht#())
DECLARE FUNCTION IClenCurt% (m%, t#(), n%, Tk#(), Ak#(), Res#)
'************************************************************************
'* 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.141592653571974                                        *
'* Error code =  0                                                      *
'*                                                                      *
'* (Exact value is -PI).                                                *
'* -------------------------------------------------------------------- *
'* REF.: "Numerical Algorithms with C, by Gisela Engeln-Muellges        *
'*        and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].           *
'*                                                                      *
'*                                Basic Release By J-P Moreau, Paris.   *
'*                                        (www.jpmoreau.fr)             *
'************************************************************************
'Program Test_Clencurt
DEFDBL A-H, O-Z
DEFINT I-N

  ISIZE = 10

  DIM t(ISIZE), Tk(ISIZE), Ak(ISIZE)

  PI = 4# * ATN(1#)

  m = 4: n = 6

  t(0) = 0#: t(1) = PI / 4#: t(2) = PI / 2#: t(3) = 3# * PI / 4#: t(4) = PI

  ierror = IClenCurtStGew(n, Tk(), Ak())

  ierror = IClenCurt(m, t(), n, Tk(), Ak(), Res)

  CLS
  PRINT
  PRINT " Integral = "; Res
  PRINT " Error code = "; ierror
  PRINT

END 'of main program

'Function to integrate
FUNCTION func (x)
  func = COS(x) - x * SIN(x)
END FUNCTION

FUNCTION IClenCurt (m, t(), n, Tk(), Ak(), Res)
'************************************************************************
'* 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:                                                           *
'*   integer m             number of subintervals                       *
'*   double  t(0:n)        partition                                    *
'*   integer n             n + 1 = number of nodes, n > 1, n even       *
'*                                   (n + 2 = global error order)       *
'*   double Tk(0:n)        Chebyshev nodes                              *
'*   double Ak(0:n)        Weights                                      *
'*                           Tk and Ak must be made available before    *
'*                           calling this function by the procedure     *
'*                           ClenCurtStGew for example                  *
'*   double Res            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                     *
'************************************************************************
  IF n < 2 OR n MOD 2 <> 0 THEN              ' n positive ? n even  ?
    IClenCurt = 1
    EXIT FUNCTION
  END IF
  IF m < 1 THEN                              ' partition
    IClenCurt = 2
    EXIT FUNCTION
  END IF

  Res = 0#
  FOR j = 0 TO m - 1                         ' loop over intervals
    v = .5# * (t(j + 1) - t(j))              ' half the interval size
    h = .5# * (t(j + 1) + t(j))              ' Interval center

    sum = 0#
    FOR k = 0 TO n                           ' Chebyshev loop
      sum = sum + Ak(k) * func(v * Tk(k) + h)
    NEXT k
    Res = Res + v * sum
  NEXT j
  IClenCurt = 0

END FUNCTION

FUNCTION IClenCurtStGew (n, StStelle(), 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:                                                          *
'*   integer n            n+1 is the number of nodes and weights:       *
'*                        n > 1, n odd                                  *
'*   double StStelle()    computed nodes                                *
'*   double Gewicht()     computed weights                              *
'*                                                                      *
'* Return value :                                                       *
'*   0:                   all ok                                        *
'*   1:                   n too small or odd                            *
'*                                                                      *
'* Author:                Uli Eggermann,  9.24.1990                     *
'************************************************************************
  DIM i AS DOUBLE

  PI = 4# * ATN(1#)

  i = 1#
  IF n < 2 OR n MOD 2 <> 0 THEN                    ' n - Test
    IClenCurtStGew = 1
    EXIT FUNCTION
  END IF
  m = n / 2
  d = 1# * (n * n - 1)
  h = 2# / (n * d)
  f = 4# / n
  p = PI / n

  Gewicht(0) = 1# / d: Gewicht(n) = 1# / d          '  left end point
  StStelle(0) = -1#: StStelle(n) = 1#               ' right end point

  FOR k = 1 TO m                                    ' Interval [-1,0]
    g = 0#
    FOR j = 1 TO m - 1
      g = g + COS(2# * j * k * p) / (4# * j * j - 1#)
    NEXT j
    Gewicht(k) = h * (d + i) - f * g
    StStelle(k) = -COS(k * p)
    i = -i
  NEXT k

  StStelle(m) = 0#                               ' center of interval

  FOR k = m + 1 TO n - 1                         ' Interval [0, 1]
    Gewicht(k) = Gewicht(n - k)
    StStelle(k) = -StStelle(n - k)
  NEXT k
  IClenCurtStGew = 0

END FUNCTION 'IClenCurtStGew