// Compile:
// g++ -std=c++0x std11_gsl_fft_arma_main.cpp `gsl-config --libs`
#include <iostream>
#include <cmath>
#include <complex>
#include <fstream>
#define ARMA
#include "gsl_fft.hpp"
#include <iomanip>
#include <armadillo>

using namespace std;
using namespace arma;

// output help routine; works with many complex data types
template <class T>
void vector_out(std::ostream& stream, const T& v){
  stream<<setprecision(8)<<scientific;
  int i=0;
  for(auto x: v){ stream<<i++<<" "<<x.real()<<" "<<x.imag()<<endl; }
}








// short name for data type
typedef Col<complex<double>> vectype;


int main (void){
  const int n=128;
  const int n1= 10; // Keep less than n!
  vectype data(n);
  ofstream myfile("result");
  //
  // Symmetric pulse:  1111111111100...001111111111
  //                   11 ones            10 ones
  data.ones();
  for(int i=n1+1;i<n-n1-1;++i) data(i)=0.0;

  cout<<"before fft :\n";
  vector_out(cout,data);
  vector_out(myfile,data);
  myfile<<endl;
  // ----------------------
  // forward FFT
  int status = fft(data,1);
  // ----------------------
  if(status!=0){
    cout << "fft fails\n";
    return 1;
  }
  cout<<"after forward fft (divided by sqrt(n)) :\n";
  data /= sqrt(n);
  vector_out(cout,data);
  vector_out(myfile,data);

  // ----------------------
  // backward FFT
  status = fft(data,-1);
  // ----------------------
  if(status!=0){
    cout << "fft fails\n";
    return 1;
  }

  cout<<"after backward fft (divided by sqrt(n)) :\n";
  data /= sqrt(n);
  vector_out(cout,data);

  myfile.close();
  return 0;
}