'********************************************************
'* Chebyshev Approximation of a user defined real       *
'* function FUNC(X) in double precision.                *
'* ---------------------------------------------------- *
'* SAMPLE RUN:                                          *
'* (Approximate sin(x) from x=0 to x=PI).               *
'*                                                      *
'* Chebyshev coefficients (N= 10 ):                     *
'*  0.944002431536470                                   *
'* -0.000000000000000                                   *
'* -0.499403258270407                                   *
'* -0.000000000000000                                   *
'*  0.027992079617546                                   *
'* -0.000000000000000                                   *
'* -0.000596695195801                                   *
'*  0.000000000000000                                   *
'*  0.000006704175524                                   *
'*  0.000000000000000                                   *
'*      X      Chebyshev Eval.   SIN(X)                 *
'*-----------------------------------------             *
'*  0.00000000   0.00000005   0.00000000                *
'*  0.34906585   0.34202018   0.34202014                *
'*  0.69813170   0.64278757   0.64278761                *
'*  1.04719755   0.86602545   0.86602540                *
'*  1.39626340   0.98480773   0.98480775                *
'*  1.74532925   0.98480773   0.98480775                *
'*  2.09439510   0.86602545   0.86602540                *
'*  2.44346095   0.64278757   0.64278761                *
'*  2.79252680   0.34202018   0.34202014                *
'*  3.14159265   0.00000005   0.00000000                *
'*                                                      *
'*                 Basic Release By J-P Moreau, Paris.  *
'*                          (www.jpmoreau.fr)           *
'* ---------------------------------------------------- *
'* REFERENCE: "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

ZERO = 0#
PI = 4# * ATN(1#)

DIM C(N), F(N)

X0 = ZERO: X1 = PI

GOSUB 1000  'call CHEBFT(X0,X1,C,N)

PRINT " Chebyshev coefficients (N="; N; "):"
F0$ = "####.###############"

FOR I = 1 TO N
 PRINT USING F0$; C(I)
NEXT I

DX = (X1 - X0) / (N - 1)
X = X0 - DX

PRINT "       X      Chebyshev Eval.   SIN(X)    "
PRINT " -----------------------------------------"

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

FOR I = 1 TO N
  X = X + DX
  GOSUB 2000  'call CHEBEV(X0,X1,C,N,X)
  PRINT USING F1$; X;
  PRINT USING F1$; CHEBEV;
  PRINT USING F1$; SIN(X)
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,C,M,X)
'**********************************************************
'* Chebyshev evaluation: All arguments are input. C is an *
'* array of Chebyshev coefficients, of length M, the first*
'* M elements of Coutput from subroutine CHEBFT (which    *
'* must have been called with the same A and B). The Che- *
'* byshev polynomial is evaluated at a point Y determined *
'* from X, A and B, and the result FUNC(X) is returned as *
'* the function value.                                    *
'**********************************************************
  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 + C(J)
    DD = SV
  NEXT J
  CHEBEV = Y * D - DD + HALF * C(1)
RETURN

'end of file Tchebysh.bas