/*********************************************************
* This program calculates the Fourier coefficients of a  *
* periodic discrete function F(x) by using the procedure *
* discreteFourierHn of unit Fourier.                     *
* ------------------------------------------------------ *
* SAMPLE RUN:                                            *
*                                                        *
* F(x) has the following shape:                          *
*                                                        *
*  1 ---- 2  pi/4                                        *
*    !  !                                                *
*    !  !                                                *
* ------!-------------------------------------------> x  *
*  -pi  !    0         ! pi                              *
*       !              !                                 *
*     3 ---------------- 4 -pi/4                         *
*                                                        *
* Calculate the Fourier coefficients of a periodic       *
* discrete function:                                     *
*                                                        *
* Number of points:  4                                   *
*                                                        *
*  x1 = -3.1415927                                       *
*  y1 = 0.7853982                                        *
*                                                        *
*  x2 = -1.5707963                                       *
*  y2 = 0.7853982                                        *
*                                                        *
*  x3 = -1.5707963                                       *
*  y3 = -0.7853982                                       *
*                                                        *
*  x4 = 3.1415927                                        *
*  y4 = -0.7853982                                       *
*                                                        *
* Lowest  harmonic: 0                                    *
* Highest harmonic: 10                                   *
*                                                        *
* a0 = -0.392699087                                      *
* b0 =  0.000000000                                      *
*                                                        *
* a1 = -0.500000023                                      *
* b1 = -0.500000048                                      *
*                                                        *
* a2 = 0.0000000025                                      *
* b2 = 0.5000000023                                      *
* -----------------                                      *
* a9 = -0.055555558                                      *
* b9 = -0.055555583                                      *
*                                                        *
* a10 = -0.000000025                                     *
* b10 =  0.100000005                                     *
*                                                        *
*                     C++ version by J-P Moreau, Paris.  *
*                             (www.jpmoreau.fr)          *
**********************************************************
Uses files: basis_r.cpp, vmblock.cpp, fourier.cpp,
            basis.h,     vmblock.h
-------------------------------------------------       */
#include <basis.h>     // ABS, MACH_EPS, REAL etc.   
#include <vmblock.h>   // dynamic allocations

int   count,h1,h2,i,npoints;
REAL  *Xi,*Yi;
REAL  ai,bi;
char  answer[3];

void  *vmblock;

//see file fourier.cpp
void discreteFourierHn(int,REAL *, REAL *, int, REAL *, REAL *);


int main() {
  
  printf("\n Calculate the Fourier coefficients of a periodic discrete function:\n\n");
  printf(" Number of points: "); scanf("%d", &npoints);
  printf("\n\n");

  // allocate vectors Xi, Yi
      vmblock = vminit();
      Xi   = (REAL *)vmalloc(vmblock, VEKTOR,  npoints+1, 0);
      Yi   = (REAL *)vmalloc(vmblock, VEKTOR,  npoints+1, 0);
      if (!vmcomplete(vmblock))
	  {
        LogError ("Memory full !", 0, __FILE__, __LINE__);
        return 1;
	  }	  

  for (i=1; i<=npoints; i++) {
    printf("  x%d = ", i); scanf("%lf", &Xi[i]);
    printf("  y%d = ", i); scanf("%lf", &Yi[i]);
    printf("\n");
  }
  printf("\n Lowest  harmonic: "); scanf("%d", &h1);
  printf("\n Highest harmonic: "); scanf("%d", &h2);
  count=1;
  for (i=h1; i<=h2; i++) {
    discreteFourierHn(npoints,Xi,Yi,i,&ai,&bi);
    printf("\n  a%d = %12.9f", i, ai);
    printf("\n  b%d = %12.9f\n", i, bi);
    if (count % 5 == 0) {
      count=0;
      scanf("%s",answer);
    }
    count++;
  }
  printf("\n\n");
  free(vmblock);
  return 0;
}

// end of file discfour.cpp