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

using namespace std;

typedef complex<double> dtype ;

void vector_out(ostream& stream, vector<dtype> v){
  vector<dtype>::iterator vpos;
  stream<<setprecision(8);
  stream<<scientific;
  int i=0;
  for(vpos=v.begin();vpos!=v.end();++vpos){
    stream<<i++<<" "<<(*vpos).real()<<" "<<(*vpos).imag()<<endl;
  }

}


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

  // pulse: 
  //     begin       nones+1                      end-nones        end
  //1    -------------------                      -----------------
  //0                       ----------------------
  transform(data.begin(),data.begin()+nones+1,data.begin(),bind2nd(plus<dtype>(),1.0));
  transform(data.end()-nones,data.end(),data.end()-nones,bind2nd(plus<dtype>(),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";
  transform(data.begin(),data.end(),data.begin(),bind2nd(divides<dtype>(),sqrt(data.size())));  
  vector_out(cout,data);
  vector_out(myfile,data);

  // ----------------------
  // backward FFT
  status = fft(data,-1);
  // ----------------------
  if(status!=0){
    cout << "fft fails\n";
    return 1;
  }
  transform(data.begin(),data.end(),data.begin(),bind2nd(divides<dtype>(),sqrt(data.size())));
  cout<<"after backward fft (divided by sqrt(n)) :\n";
  vector_out(cout,data);

  myfile.close();
  return 0;
}