#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <fftw3.h>
#include <cmath>

#include "JLang/gzstream.h"

#include "Jeep/JeepToolkit.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 * Auxiliary program to run fast Fourier transformation (FFT).\n
 * For more information, see <a href="https://www.fftw.org/">FFTW</a>.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string  inputFile;
  string  outputFile;
  double  binWidth;
  int     precision;
  int     debug;

  try { 

    JParser<> zap("Auxiliary program to run fast Fourier transformation (FFT).");
    
    zap['f'] = make_field(inputFile,  "input file (containing 1D array of values)");
    zap['o'] = make_field(outputFile, "output file (containing 2D array of values)")   = "fft.txt";
    zap['B'] = make_field(binWidth,   "bin width of input array (if zero, use index)") = 0.0;
    zap['F'] = make_field(precision,  "number of decimals if output format")           = 9;
    zap['d'] = make_field(debug)           = 2;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }

  // input
  
  vector<double> data;

  istream* in = open<istream>(inputFile);

  if (in == NULL) {
    FATAL("Invalid file " << inputFile << endl);
  }

  for (double x; *in >> x; ) {
    data.push_back(x);
  }

  close(in);

  const size_t N  = data.size(); 
 
  if (N == 0) {
    FATAL("Number of points: " << N << endl);
  }
  
  // FFT

  const size_t  NS     = N/2 + 1;
  double*       input  = (double*) malloc(sizeof(double) *  N);
  fftw_complex* output = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * NS);
  fftw_plan     plan   = fftw_plan_dft_r2c_1d(N, input, output, FFTW_MEASURE);

  for (size_t i = 0; i != N; ++i) {
    input[i] = data[i];
  }

  // execute plan

  fftw_execute(plan);
	
  // destroy plan

  fftw_destroy_plan(plan);

  // output

  const double xmin  = (binWidth != 0.0 ? 1.0 / (N * binWidth) : 0.0);
  const double xmax  = (binWidth != 0.0 ? 1.0 / (binWidth)     : 0.0);

  const int    width = precision + 7;

  ostream* out = open<ostream>(outputFile.c_str());

  for (size_t i = 0; i != NS; ++i) {

    const double real  = output[i][0];
    const double imag  = output[i][1];
    const double power = sqrt(real*real + imag*imag);

    if (binWidth != 0.0)
      *out << SCIENTIFIC(width,precision) << xmin + (double) i * (xmax - xmin) / (double) (N - 1);
    else
      *out << i;

    *out << ' ';
    *out << SCIENTIFIC(width,precision) << power << endl;
  }

  fftw_free(output);
  free(input);

  close(out);  
}