#include <math.h>
#include <iostream.h>

#include "fftpserver.h"
#include "sambainit.h"

#include "timing.h"


template <class T>
inline T module(complex<T> c) 
{
  return (sqrt(c.real()*c.real()+c.imag()*c.imag()));
}

template <class T>
void TestFFTS(T seuil, FFTServerInterface & ffts)
{
  
  int num=16;
  int i;

  T fact = 1./num;

  TVector< complex<T> > inc(num), bkc(num), difc(num);
  TVector< T > in(num), ino(num), bk(num),dif(num);
  TVector< complex<T> > outc(num);
  
  
  for (int i=0; i<num ; i++){
    ino[i] = in[i] = cos(2*M_PI*i/8.) + sin(2*M_PI*i/4.);
    inc[i] = complex<T> (in[i], 0.);
  }
  

  cout << " Testing FFTServer " << ffts.getInfo() << endl;

  cout << "Input / L = " << num << endl; 
  in.Print(0,0,num,0,num);
  cout << endl;

  cout << " >>>> Testing FFTPackServer "  << endl;
  FFTPackServer fftp;
  cout << " Testing FFTPackServer "  << endl;
  fftp.fftf(in.NElts(), in.Data());
  in /= (num/2.);
  cout << " fftp.fftf(in.NElts(), in.Data()) FORWARD: " << endl; 
  in.Print(0,0,num,0,num);
  cout << endl;
  fftp.fftb(in.NElts(), in.Data());
  cout << " fftp.fftb(in.NElts(), in.Data()) BACKWARD: " << endl; 
  in.Print(0,0,num,0,num);
  cout << endl;

  in = ino;

  cout << "\n ----  Testing FFT(T, complex<T>) ---- " << endl;

  ffts.FFTForward(in, outc);
  cout << " outc , NElts= " << outc.NElts() << endl;
  outc.Print(0,0,num,0,num);

  ffts.FFTBackward(outc, bk);
  cout << " bk , NElts= " << bk.NElts() << endl;
  bk.Print(0,0,num,0,num);
 
  dif = bk*fact - in;
  cout << " dif , NElts= " << dif.NElts() << endl;
  dif.Print(0,0,num,0,num);

  int ndiff = 0;
  for(i=0; i<num; i++) 
    if (fabs(dif(i)) > seuil)  ndiff++;

  cout << " Difference, Seuil= " << seuil << " NDiff= " << ndiff << endl;
 
  cout << "\n ----  Testing FFT(complex<T>, complex<T>) ---- " << endl;
  ffts.FFTForward(inc, outc);
  cout << " outc , NElts= " << outc.NElts() << endl;
  outc.Print(0,0,num,0,num);

  ffts.FFTBackward(outc, bkc);
  cout << " bkc , NElts= " << bkc.NElts() << endl;
  bkc.Print(0,0,num,0,num);
 
  difc = bkc*complex<T>(fact,0.) - inc;
  cout << " difc , NElts= " << difc.NElts() << endl;
  difc.Print(0,0,num,0,num);

  int ndiffc = 0;
  for(i=0; i<num; i++) 
    if (fabs(module(difc(i))) > seuil)  ndiffc++;
  cout << " Difference, Seuil= " << seuil << " NDiffC= " << ndiffc << endl;
}


int main(int narg, char* arg)
{

  PeidaInit();
  InitTim();   // Initializing the CPU timer
  FFTPackServer ffts;


  try {
  float fs = 1.e-6;
  double ds = 1.e-12;
  cout << "\n ======================================== \n"
       << "    ------ Testing FFTS for float ----- \n"
       << " ======================================== " << endl;
  TestFFTS(fs, ffts);
  cout << "\n ======================================== \n"
       << "    ------ Testing FFTS for double ----- \n"
       << " ======================================== " << endl;
  TestFFTS(ds, ffts);
  }
  catch(PThrowable exc ) {
    cerr << "TestFFT-main() , Catched exception: " << exc.Msg() << endl;
  }
  catch(std::exception ex) {
    cerr << "TestFFT-main() , Catched exception ! " << (string)(ex.what()) << endl;
  }
  /*
  catch(...) {
    cerr << "TestFFT-main() , Catched ... exception ! " << endl;
  }
  */
  PrtTim("End of tfft ");
}