'******************************************************
'* Chebyshev Approximation of a user defined real     *
'* function FUNC(X) in double precision.              *
'* -------------------------------------------------- *
'* SAMPLE RUN:                                        *
'* (Approximate integral of sin(x) from x=0 to x=PI). *
'*                                                    *
'*      X          Chebyshev Eval.   -COS(X)+1        *
'*                  of integral                       *
'*------------------------------------------------    *
'*  0.00000000      -0.00000000      0.00000000       *
'*  0.15707963       0.01231166      0.01231166       *
'*  0.31415927       0.04894348      0.04894348       *
'*  0.47123890       0.10899348      0.10899348       *
'*  0.62831853       0.19098301      0.19098301       *
'*  0.78539816       0.29289321      0.29289322       *
'*  0.94247780       0.41221474      0.41221475       *
'*  1.09955743       0.54600950      0.54600950       *
'*  1.25663706       0.69098301      0.69098301       *
'*  1.41371669       0.84356554      0.84356553       *
'*  1.57079633       1.00000000      1.00000000       *
'*  1.72787596       1.15643446      1.15643447       *
'*  1.88495559       1.30901699      1.30901699       *
'*  2.04203522       1.45399050      1.45399050       *
'*  2.19911486       1.58778526      1.58778525       *
'*  2.35619449       1.70710679      1.70710678       *
'*  2.51327412       1.80901699      1.80901699       *
'*  2.67035376       1.89100652      1.89100652       *
'*  2.82743339       1.95105651      1.95105652       *
'*  2.98451302       1.98768834      1.98768834       *
'*  3.14159265       2.00000000      2.00000000       *
'*                                                    *
'*                Basic Release By J-P Moreau, Paris. *
'*                         (www.jpmoreau.fr)          *
'* -------------------------------------------------- *
'* REF.:  "Numerical Recipes, The Art of Scientific   *
'*         Computing By W.H. Press, B.P. Flannery,    *
'*         S.A. Teukolsky and W.T. Vetterling,        *
'*         Cambridge University Press, 1986"          *
'*         [BIBLI 08].                                *
'******************************************************
'PROGRAM TESTCHEBY
DEFDBL A-H, O-Z
DEFINT I-N

N = 10

ONE = 1#
QUART = .25#
ZERO = 0#

PI = 4# * ATN(1#)

DIM C(N), CIN(N), F(N)

X0 = ZERO: X1 = PI

'calculate Chebyshev coefficients of FUNC(X)
GOSUB 1000  'call CHEBFT(X0,X1,C,N)

'calculate Chebyshev coefficients of integral of FUNC(X)
GOSUB 3000  'call CHINT(X0,X1,C,CIN,N)

NPTS = 2 * N + 1
DX = (X1 - X0) / (NPTS - 1)
X = X0 - DX

PRINT "       X      Chebyshev Eval.  -COS(X)+1  "
PRINT "               of integral                "
PRINT " -----------------------------------------"

F1$ = "####.########"

FOR I = 1 TO NPTS
  X = X + DX
  GOSUB 2000  'call CHEBEV(X0,X1,CIN,N,X)
  PRINT USING F1$; X;
  PRINT USING F1$; CHEBEV;
  PRINT USING F1$; -COS(X) + ONE
NEXT I

END 'of main program

'user defined function FUNC(XX)
500 FUNC = SIN(XX)
RETURN

1000 'subroutine CHEBFT(A,B,C,N)
'********************************************************
'* Chebyshev fit: Given a real function FUNC(X), lower  *
'* and upper limits of the interval [A,B] for X, and a  *
'* maximum degree N, this routine computes the N Cheby- *
'* shev coefficients Ck, such that FUNC(X) is approxima-*
'* ted by:  N                                           *
'*         [Sum Ck Tk-1(Y)] - C1/2, where X and Y are   *
'*         k=1                                          *
'* related by:     Y = (X - 1/2(A+B)) / (1/2(B-A))      *
'* This routine is to be used with moderately large N   *
'* (e.g. 30 or 50), the array of C's subsequently to be *
'* truncated at the smaller value m such that Cm+1 and  *
'* subsequent elements are negligible.                  *
'********************************************************
  HALF = .5#: TWO = 2#
  A = X0: B = X1
  BMA = HALF * (B - A): BPA = HALF * (B + A)
  FOR K = 1 TO N
    Y = COS(PI * (K - HALF) / N)
    XX = Y * BMA + BPA: GOSUB 500
    F(K) = FUNC
  NEXT K
  FAC = TWO / N
  FOR J = 1 TO N
    SUM = ZERO
    FOR K = 1 TO N
      SUM = SUM + F(K) * COS((PI * (J - 1)) * ((K - HALF) / N))
    NEXT K
    C(J) = FAC * SUM
  NEXT J
RETURN


2000 'function CHEBEV(A,B,CIN,M,X)
'*********************************************************
'* Chebyshev evaluation: All arguments are input. CIN is *
'* an array of Chebyshev coefficients, of length M, the  *
'* first M elements of C output from subroutine CHINT    *
'* (which must have been called with the same A and B).  *
'* The Chebyshev polynomial is evaluated at a point Y    *
'* determined from X, A and B, and the result, integral  *
'* of FUNC(X), is returned in variable CHEBEV.           *
'*********************************************************
  IF (X - A) * (X - B) > ZERO THEN
    PRINT " X not in range."
    RETURN
  END IF
  M = N
  D = ZERO: DD = ZERO
  Y = (TWO * X - A - B) / (B - A)'change of variable
  Y2 = TWO * Y
  FOR J = M TO 2 STEP -1
    SV = D
    D = Y2 * D - DD + CIN(J)
    DD = SV
  NEXT J
  CHEBEV = Y * D - DD + HALF * CIN(1)
RETURN

3000 'subroutine CHINT(A,B,C,CIN,N)
'**********************************************************
'* Given A,B,C, as output from routine CHEBFT, and given  *
'* N, the desired degree of approximation (length of C to *
'* be used), this routine returns the array CIN, the Che- *
'* byshev coefficients of the integral of the function    *
'* whose coefficients are C. The constant of integration  *
'* is set so that the integral vanishes at A.             *
'**********************************************************
  CON = QUART * (B - A)
  SUM = ZERO
  FAC = ONE
  FOR J = 2 TO N - 1
    CIN(J) = CON * (C(J - 1) - C(J + 1)) / (J - 1)
    SUM = SUM + FAC * CIN(J)
    FAC = -FAC
  NEXT J
  CIN(N) = CON * C(N - 1) / (N - 1)
  SUM = SUM + FAC * CIN(N)
  CIN(1) = TWO * SUM      'set the constant of integration
RETURN

'end of file Tchebint.bas