'************************************************************************
'*                                                                      *
'*           Interpolate a function by trigonometric polynoms           *
'*                                                                      *
'* -------------------------------------------------------------------- *
'*   SAMPLE RUN:                                                        *
'*                                                                      *
'*   Number of points: 3                                                *
'*                                                                      *
'*   Number of harmonics: 3                                             *
'*                                                                      *
'*                                                                      *
'*   Function to interpolate:                                           *
'*      X       Y                                                       *
'*     0.000000   3.141593                                              *
'*     0.897598   0.897598                                              *
'*     1.795196   1.795196                                              *
'*     2.692794   2.692794                                              *
'*     3.590392   3.590392                                              *
'*     4.487990   4.487990                                              *
'*     5.385587   5.385587                                              *
'*     6.283185   3.141593                                              *
'*                                                                      *
'*   Fourier coefficients:                                              *
'*   A( 0 ) =    3.141593  B( 0 ) =    0.000000                         *
'*   A( 1 ) =   -0.000000  B( 1 ) =   -1.863881                         *
'*   A( 2 ) =   -0.000000  B( 2 ) =   -0.715810                         *
'*   A( 3 ) =   -0.000000  B( 3 ) =   -0.204871                         *
'*                                                                      *
'*   Interpolation:                                                     *
'*      X       Y                                                       *
'*     0.000000   3.141593                                              *
'*     0.448799   1.573507                                              *
'*     0.897598   0.897598                                              *
'*     1.346397   1.174039                                              *
'*     1.795196   1.795196                                              *
'*     2.243995   2.293325                                              *
'*     2.692794   2.692794                                              *
'*     3.141593   3.141593                                              *
'*     3.590392   3.590392                                              *
'*     4.039191   3.989860                                              *
'*     4.487990   4.487990                                              *
'*     4.936788   5.109147                                              *
'*     5.385587   5.385587                                              *
'*     5.834386   4.709678                                              *
'*     6.283185   3.141593                                              *
'*                                                                      *
'* -------------------------------------------------------------------- *
'*    "Méthodes de calcul numérique, Tome 2 By Claude Nowakowski,       *
'*     PSI Edition, 1984" [BIBLI 04].                                   *
'*                                                                      *
'*                             Basic Release 1.0 By J-P Moreau, Paris.  *
'*                                       (www.jpmoreau.fr)              *
'************************************************************************
'Program TrigPoly
DEFDBL A-H, O-Z
DEFINT I-N

PI = 3.1415926535#

CLS
10 PRINT
INPUT " Number of points: ", N
PRINT
INPUT " Number of harmonics: ", NH
PRINT
IF N < NH THEN
  PRINT " Error, N must be >= NH."
  GOTO 10
END IF

NN = 2 * N + 1

' TF stores function values
DIM TF(NN)
FOR K = 0 TO NN
  TF(K) = K * 2# * PI / NN
NEXT K

TF(0) = (TF(0) + TF(NN)) / 2#
TF(NN) = TF(0)

CF = 2# / NN

F$ = "####.######"

PRINT
PRINT " Function to interpolate:"
PRINT "    X       Y"
FOR K = 0 TO NN
  X = K * 2# * PI / NN
  PRINT USING F$; X;
  PRINT USING F$; TF(K)
NEXT K

CT = PI * CF

'Initialize
S1 = SIN(CT): C1 = COS(CT)
CP = 1#: SP = 0#: IP = 0: F0 = TF(0)

DIM A(NH + 1), B(NH + 1)'to store Fourier coeff.

'Calculate Fourier coefficients A, B, for each harmonic
90  U2 = 0#: U1 = 0#: M = NN - 1
100 U0 = TF(M) + 2# * CP * U1 - U2
U2 = U1: U1 = U0: M = M - 1
IF M > 0 THEN GOTO 100

A(IP) = CF * (F0 + CP * U1 - U2)
B(IP) = CF * SP * U1

IF IP > NH THEN GOTO 200
Q = C1 * CP - S1 * SP
SP = C1 * SP + S1 * CP
CP = Q
IP = IP + 1
GOTO 90

200 A(0) = .5# * A(0)
PRINT
PRINT " Fourier coefficients:"
FOR K = 0 TO NH
  PRINT " A("; K; ") = ";
  PRINT USING F$; A(K);
  PRINT "  B("; K; ") = ";
  PRINT USING F$; B(K)
NEXT K

'Interpolate function
PRINT
PRINT " Interpolation:"
PRINT "    X       Y"
FOR K = 0 TO NN * 2
  X = K * 2# * PI / (2# * NN)
  S = A(0)
  FOR L = 1 TO NH
    S = S + A(L) * COS(X * L) + B(L) * SIN(X * L)
  NEXT L
  PRINT USING F$; X;
  PRINT USING F$; S
NEXT K

END 'of main program