'******************************************************
'* Chebyshev Approximation of a user defined real     *
'* function FUNC(X) in double precision.              *
'* -------------------------------------------------- *
'* SAMPLE RUN:                                        *
'* (Approximate derivative of sin(x) from x=0 to      *
'*  x=PI).                                            *
'*                                                    *
'*      X          Chebyshev Eval.    COS(X)          *
'*                 of derivative                      *
'*-----------------------------------------------     *
'*  0.00000000      0.99999707      1.00000000        *
'*  0.15707963      0.98768900      0.98768834        *
'*  0.31415927      0.95105644      0.95105652        *
'*  0.47123890      0.89100611      0.89100652        *
'*  0.62831853      0.80901694      0.80901699        *
'*  0.78539816      0.70710708      0.70710678        *
'*  0.94247780      0.58778552      0.58778525        *
'*  1.09955743      0.45399047      0.45399050        *
'*  1.25663706      0.30901672      0.30901699        *
'*  1.41371669      0.15643421      0.15643447        *
'*  1.57079633      0.00000000      0.00000000        *
'*  1.72787596     -0.15643421     -0.15643447        *
'*  1.88495559     -0.30901672     -0.30901699        *
'*  2.04203522     -0.45399047     -0.45399050        *
'*  2.19911486     -0.58778552     -0.58778525        *
'*  2.35619449     -0.70710708     -0.70710678        *
'*  2.51327412     -0.80901694     -0.80901699        *
'*  2.67035376     -0.89100611     -0.89100652        *
'*  2.82743339     -0.95105644     -0.95105652        *
'*  2.98451302     -0.98768900     -0.98768834        *
'*  3.14159265     -0.99999707     -1.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 TESTCHDER
DEFDBL A-H, O-Z
DEFINT I-N

N = 10

QUART = .25#
ZERO = 0#

PI = 4# * ATN(1#)

DIM C(N), CDER(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 derivative of FUNC(X)
GOSUB 3000  'call CHDER(X0,X1,C,CDER,N)

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

CLS
PRINT
PRINT "       X       Chebyshev Eval.  COS(X)   "
PRINT "               of derivative             "
PRINT " ----------------------------------------"

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

FOR I = 1 TO NPTS
  X = X + DX
  GOSUB 2000  'call CHEBEV(X0,X1,CDER,N,X)
  PRINT USING F1$; X;
  PRINT USING F1$; CHEBEV;
  PRINT USING F1$; COS(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,CDER,M,X)
'*********************************************************
'* Chebyshev evaluation: All arguments are input. CDER is*
'* an array of Chebyshev coefficients, of length M, the  *
'* first M elements of CDER output from subroutine 3000  *
'* (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 + CDER(J)
    DD = SV
  NEXT J
  CHEBEV = Y * D - DD + HALF * CDER(1)
RETURN

3000 'subroutine CHDER(A,B,C,CDER,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 CDER, the Che-*
'* byshev coefficients of the derivative of the function  *
'* whose coefficients are C.                              *
'**********************************************************
CDER(N) = ZERO
CDER(N - 1) = TWO * (N - 1) * C(N)
IF N >= 3 THEN
  FOR J = N - 2 TO 1 STEP -1
    CDER(J) = CDER(J + 2) + TWO * J * C(J + 1)
  NEXT J
END IF
CON = TWO / (B - A)
FOR J = 1 TO N            'normalize to interval B - A
  CDER(J) = CDER(J) * CON
NEXT J
RETURN

'end of file Tchebder.bas