/***********************************************************************
*                                                                      *
*           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].                                   *
*                                                                      *
*                               C++ Release 1.0 By J-P Moreau, Paris.  *
*                                         (www.jpmoreau.fr)            *
***********************************************************************/
#include <stdio.h>
#include <math.h>

#define  SIZE  50
#define  PI    3.1415926535

// Labels: e10, e90, e100, e200

    int IP, K, L, M, N, NH, NN;
    double A[SIZE], B[SIZE], TF[SIZE];
    double C1, CF, CP, CT, F0, S1, SP, X;
    double Q, S, U0, U1, U2;


void main() {

e10:printf("\n Number of points   : "); scanf("%d", &N);
    printf("\n Number of harmonics: "); scanf("%d", &NH);
    printf("\n");

    if (N < NH) {
      printf(" Error, N must be >= NH.\n");
      goto e10;
    }

    NN = 2 * N + 1;

// TF stores function values
    for (K=0; K<=NN; K++)
      TF[K] = K * 2.0 * PI / NN;

    TF[0]  = (TF[0] + TF[NN]) / 2.0;
    TF[NN] = TF[0];

    CF = 2.0 / NN;

    printf("\n Function to interpolate:\n");
    printf("      X          Y\n");
    for (K=0; K<=NN; K++) {
      X = K * 2.0 * PI / NN;
      printf("%11.6f%11.6f\n", X, TF[K]);
    }

    CT = PI * CF;

// Initialize
    S1 = sin(CT); C1 = cos(CT);
    CP = 1.0; SP = 0.0; IP = 0; F0 = TF[0];

// Calculate Fourier coefficients A, B, for each harmonic
e90: U2 = 0.0; U1 = 0.0; M = NN - 1;
e100:U0 = TF[M] + 2.0 * CP * U1 - U2;
    U2 = U1; U1 = U0; M--;
    if (M > 0) goto e100;

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

    if (IP > NH) goto e200;
    Q = C1 * CP - S1 * SP;
    SP = C1 * SP + S1 * CP;
    CP = Q;
    IP++;
    goto e90;

e200:A[0] = 0.5 * A[0];
    printf("\n Fourier coefficients:\n");
    for (K=0; K<=NH; K++)
      printf(" A(%d) = %11.6f  B(%d) = %11.6f\n", K, A[K], K, B[K]);

// Interpolate function
    printf("\n Interpolation:\n");
    printf("      X          Y\n");
    for (K=0; K<=NN*2; K++) {
      X = K * 2.0 * PI / (2 * NN);
      S = A[0];
      for (L=1; L<=NH; L++)
        S += A[L] * cos(X * L) + B[L] * sin(X * L);
      printf("%11.6f%11.6f\n", X, S);
    }

    printf("\n");

} // end of main program

// End of file trigpoly.cpp