/* gsl_inv_fft.cc

calculates the inverse fft using gsl fft routines.  assumes that the
fft itself was originally calculated from purely real data (The means
the fft is stored in the "half-complex" scheme, see gsl
documentation.)



related program(s): gsl_fft.cc
                    calculates the fft of a purely real function.




*/

#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 data[i_max];
       double x;
       
       gsl_fft_halfcomplex_wavetable * hc;
       gsl_fft_real_workspace * work;
       
       ifstream input_file;
       ofstream output_file;

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

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

       i = 0;

       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 actually read



       work = gsl_fft_real_workspace_alloc (N);
     
       hc = gsl_fft_halfcomplex_wavetable_alloc (N);
     
       gsl_fft_halfcomplex_inverse (data, 1, N, 
				hc, work);
       gsl_fft_halfcomplex_wavetable_free (hc);
     
       for (int i = 0;  i < N; ++i)
     {
       output_file << i << "\t" << data[i] << endl;
     }

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