// Fourier series routines
#include <math.h>

#define  MAXITER  15

double om;
double PI = 3.1415926535;
int flag;

// analytical function to study:
// ----------------------------
// Note: Exact Fourier coefficients for this periodic function are:
//  an=0 if n is even, -1/2n if n=4p+1, 1/2n if n=4p+3
//  bn=0 if n=4p, -1/2n if n=4p+1, 1/n if n=4p+2, -1/2n if n=4p+3
//----------------------------------------------------------------
  double F(double x) {
    if (x<-PI) return 0.0;
    else if (x<-PI/2)  return PI/4;
    else if (x<=PI) return -PI/4;
    else return 0.0;
  }

// Function to integrate by Romberg method
  double FUNC(double x) {
    if (!flag)  return (F(x)*cos(om*x));  //for an
    else  return (F(x)*sin(om*x));        //for bn
  }

  /***************************************************************
  *   Calculate the Fourier harmonic #n of a periodic discreet   *
  *   function F(x) defined by ndata points.                     *
  * ------------------------------------------------------------ *
  * Inputs:                                                      *
  *            ndata: number of points of discreet function.     *
  *            X    : pointer to table storing xi abscissas.     *
  *            Y    : pointer to table storing yi ordinates.     *
  *                                                              *
  * Outputs:                                                     *
  *            a    : coefficient an of the Fourier series.      *
  *            b:   : coefficient bn of the Fourier series.      *
  * ------------------------------------------------------------ *
  * Reference: "Mathematiques en Turbo-Pascal By Marc Ducamp and *
  *             Alain Reverchon, 1. Analyse, Editions Eyrolles,  *
  *             Paris, 1991" [BIBLI 03].                         *
  *                                                              *
  *                           C++ version by J-P Moreau, Paris.  *
  ****************************************************************
  Note: The Fourier series  of a periodic discreet function F(x)
        can be written under the form:
                    n=inf.
        F(x) = a0 + Sum ( an cos(2 n pi/T x) + bn sin(2 n pi/T x)  
                    n=1 
  ----------------------------------------------------------------*/  
void DiscreetFourierHn(int ndata, double *X, double *Y, int n, double *a, double *b)  {
  int i;
  double xi,T,om,wa,wb,wc,wd,wg,wh,wi,wl,wm,wn,wp;
  T=X[ndata]-X[1]; xi=(X[ndata]+X[1])/2;
  om = 2*PI*n/T; *a=0; *b=0;
  for (i=1; i<ndata; i++) {
    wa=X[i]; wb=X[i+1];
    wc=Y[i]; wd=Y[i+1];
    if (wa != wb)  {
      wg = (wd-wc)/(wb-wa);
      wh = om*(wa-xi); wi=om*(wb-xi);
      if (n==0)
        *a += (wb-wa)*(wc+wg/2*(wb-wa));
      else {
        wl = cos(wh); wm = sin(wh);
        wn = cos(wi); wp = sin(wi);
        *a += wg/om*(wn-wl) + wd*wp - wc*wm;
        *b += wg/om*(wp-wm) - wd*wn + wc*wl;
      }
    }
  }
  *a /= T; *b /= T;
  if (n != 0)  {
    *a *= 2/om; *b *= 2/om;
  }
}

  /******************************************************
  * Integral of a function FUNC(X) by Romberg's method  *
  * --------------------------------------------------- *
  * INPUTS:                                             *
  *          a      begin value of x variable           *
  *          b      end value of x variable             *
  *       prec      desired precision                   *
  *    itermin      minimum number of ietrations        *
  *    itermax      maximum number of iterations        *
  *                                                     *
  * OUTPUTS:                                            *
  *    obtprec      obtained precision for integral     *
  *          n      number of iterations done           *
  *                                                     *
  * RETURNED VALUE  the integral of FUNC(X) from a to b *
  *                                                     *
  ******************************************************/ 
  double RombergIntegral(double a,double b,double prec, double *obtprec,
                                           int *n, int itermin, int itermax)  {
    int    i,j;
    double pas,r,s,ta;
    double t[MAXITER+1][MAXITER+1];
    if (itermax>MAXITER) itermax=MAXITER;
    r = FUNC(a);
    ta = (r + FUNC(b)) / 2;
    *n=0;
    pas=b-a;
    t[0][0]=ta*pas;
	do {
      *n = *n + 1;
      pas = pas/2;
      s=ta;
	  for (i=1; i<pow(2,*n); i++)
        s += FUNC(a+pas*i);
      t[0][*n]=s*pas;
      r=1;
	  for (i=1; i<*n+1; i++) {
        r=r*4;
        j=*n-i;
        t[i][j]=(r*t[i-1][j+1] - t[i-1][j])/(r-1);
	  }
      *obtprec = fabs(t[*n][0] - t[*n-1][0]);
      if (*n>=itermax) return t[*n][0];
	} while (*obtprec >prec || *n<itermin); 
    return t[*n][0];
  }

  /***************************************************************
  * Calculate the Fourier harmonic #n of a periodic function F(x)*
  * analytically defined.                                        *
  * ------------------------------------------------------------ *
  * Inputs:                                                      *
  *            t1   : -period/2                                  *
  *            t2   : period/2                                   *
  *            n    : order of harmonic                          *          
  *                                                              *
  * Outputs:                                                     *
  *            a    : coefficient an of the Fourier series.      *
  *            b:   : coefficient bn of the Fourier series.      *
  * ------------------------------------------------------------ *
  * Reference: "Mathematiques en Turbo-Pascal by Marc Ducamp and *
  *             Alain Reverchon, 1. Analyse, Editions Eyrolles,  *
  *             Paris, 1991" [BIBLI 03].                         *
  *                                                              *
  *                           C++ version by J-P Moreau, Paris.  *
  ****************************************************************
  Note:  When a periodic function f(x), of period T, can be developped
         in a Fourier series, this one is unique:
                     n=inf.
         F(x) = a0 + Sum ( an cos(2 n pi/T x) + bn sin(2 n pi/T x)  
                     n=1
         (b0 is always zero).

         The coefficients an, bn are given by:
                      T/2
             a0 = 1/T Sum (f(x) dx)
                      -T/2
                      T/2
             an = 2/T Sum (f(x) cos(2npi/T x) dx)
                      -T/2
                      T/2
             an = 2/T Sum (f(x) sin(2npi/T x) dx)
                      -T/2
         Here, the integrals are calculated by the Romberg method.
  --------------------------------------------------------------------*/
void AnalyticFourierHn(double t1,double t2, int n, double *a, double *b)  {
  double T,precision,temp;
  int itmax,niter;
  T=fabs(t2-t1);
  om=2.0*PI*n/T;
  precision=1e-10; itmax=15;
  flag=0; // calculate an
  *a=RombergIntegral(t1,t2,precision,&temp,&niter,5,itmax)*2/T;
  flag=1; // calculate bn
  *b=RombergIntegral(t1,t2,precision,&temp,&niter,5,itmax)*2/T;
  if (n==0) *a=*a/2;
}

// end of file fourier.cpp