/*

 gsl_fft_tec.cc

code to compute the fft of  a real function
using gsl routines.

related program(s): gsl_inv_fft.cc 

*/


#include <iostream>
#include <cmath>
#include <fstream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>
     

using namespace std;

const int i_max = 1000;

int main ()
{
  int i = 0, j;
  double x;
  double data[i_max];
       
  gsl_fft_real_wavetable * real;
  //       gsl_fft_halfcomplex_wavetable * hc;
  gsl_fft_real_workspace * work;

  ifstream input_file;
  ofstream output_file;

  input_file.open("fourier_raw.dat");   // get input data
  output_file.open("gsl_fft.dat");          // save processed data


  for (i = 0; i < i_max; ++i) // zero out working array
    {
      data[i] = 0.0;
    }

  i = 0; // reset i

  while (!input_file.eof() && (i < i_max)){  // fill sampled data array
    input_file >> j >> x;
    data[i] = x;
    ++i;
  }

       int N = i;   // number of data points read

       
       work = gsl_fft_real_workspace_alloc (i_max);
       real = gsl_fft_real_wavetable_alloc (i_max);
       
       gsl_fft_real_transform (data, 1, i_max, real, work);
       
       gsl_fft_real_wavetable_free (real);
       
       for (int i = 0;  i < N; ++i)
	 {
	   output_file << i << "\t" << data[i] << endl;
	 }

       output_file.close();
       cout << " data saved in gsl_fft.dat" << endl;
       gsl_fft_real_workspace_free (work);
       return 0;
}