#include <iostream>
#include <cmath>
#include <vector>
#include <complex>
#include <iterator>
#include <algorithm>
#include <fstream>
#include "gsl_fft.hpp"
#include <iomanip>
using namespace std;

template <class T>
void vector_out(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 complex<double> dtype ;

int main (void){
  const int n=128;
  const int n1= 10; // must be less than n!
  vector<dtype> data(n,0.0);
  ofstream myfile("result");

  // pulse: 
  //     begin       n1+1                      end-n1        end
  //1    -------------------                      -----------------
  //0                       ----------------------
  for_each(data.begin(),data.begin()+n1+1,[&](dtype & d){d=1.0;}); //lambda function
  for_each(data.end()-n1,data.end(),[&](dtype & d){d=1.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";
  for( auto & x: data) {x/=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;
  }
  for( auto & x: data) {x/=sqrt(n);}
  cout<<"after backward fft (divided by sqrt(n)) :\n";
  vector_out(cout,data);

  myfile.close();
  return 0;
}