!************************************************************************
!*                                                                      *
!*           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].                                   *
!*                                                                      *
!*                               F90 Release 1.0 By J-P Moreau, Paris.  *
!*                                        (www.jpmoreau.fr)             *
!************************************************************************
Program TrigPoly

real*8, parameter :: PI = 3.1415926535d0
integer,parameter :: SIZE = 25

real*8 A(0:SIZE), B(0:SIZE), TF(0:SIZE)

real*8 C1, CF, CP, CT, F0, S1, SP, X
real*8 Q, S, U0, U1, U2

integer IP, K, L, M, N, NH, NN


10 PRINT *,' '
write(*,11,advance='no'); read *, N
PRINT *,' '
write(*,12,advance='no'); read *, NH
PRINT *,' '
IF (N < NH) THEN
  PRINT *,' Error, N must be >= NH.'
  GOTO 10
END IF

NN = 2 * N + 1

print *, N, NH, NN

! TF stores function values
DO K=0, NN
  TF(K) = K * 2.D0 * PI / DFLOAT(NN)
END DO

TF(0) = (TF(0) + TF(NN)) / 2.D0
TF(NN) = TF(0)

CF = 2.D0 / DFLOAT(NN)

PRINT *,' '
PRINT *,' Function to interpolate:'
PRINT *,'      X         Y'
DO K=0, NN
  X = K * 2.D0 * PI / DFLOAT(NN)
  write(*,20)  X, TF(K)
END DO

CT = PI * CF

!Initialize
S1 = DSIN(CT); C1 = DCOS(CT)
CP = 1.D0; SP = 0.D0; IP = 0; F0 = TF(0)

! Calculate Fourier coefficients A, B, for each harmonic
90  U2 = 0.D0; U1 = 0.D0; M = NN - 1
100 U0 = TF(M) + 2.D0 * CP * U1 - U2
U2 = U1; U1 = U0; M = M - 1
IF (M > 0) GOTO 100

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

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

200 A(0) = .5D0 * A(0)
PRINT *,' '
PRINT *,' Fourier coefficients:'
DO K=0, NH
  write(*,30) K, A(K), K, B(K)
END DO

! Interpolate function
PRINT *,' '
PRINT *,' Interpolation:'
PRINT *,'      X         Y'
DO K=0, NN*2
  X = K * 2.D0 * PI / (2.D0 * NN)
  S = A(0)
  DO L=1, NH
    S = S + A(L) * DCOS(X * L) + B(L) * DSIN(X * L)
  END DO
  write(*,20) X, S
END DO

stop

11 format(' Number of points   : ')
12 format(' Number of harmonics: ')

20 format(2F11.6)
30 format(' A(',I2,') = ',F11.6,'  B(',I2,') = ',F11.6)

END !of main program

!end of file trigpoly.f90