/*----------------------------------------------------------------
  This program show the use of 1d complex Fourier transformation using 
  FFTW3.
  
  Note that the input data must be in the following form :
  1d : x, Re(f(x)),Im(f(x))
  
  For comments and suggestions mail to prasad.jayanti@gmail.com
  Jayanti Prasad 
  Dec 04, 2013  
  -----------------------------------------------------------------*/


#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<fftw3.h> 
// fftw3.h must be included
#include<complex.h>
// complex.h is needed for complex numbers 

int modes(int , float []);

int main(int argc, char *argv[]) {
  int i, n;
  fftw_complex *X,*X1,*Y; 
  fftw_plan p,q; 
  FILE *fp1,*fp2;
  float *x,*f,t2,t3,C,f_nyq,dx; 
  
  if (argc < 2){
    fprintf(stderr,"./example_dft_1d  <input file> <size along one dim>\n");
    return(-1);
  }
  
  fp1 = fopen(argv[1],"r");
  n   = atoi(argv[2]);
  C   = 1.0/n;
  
  X  = fftw_alloc_complex(n);
  X1 = fftw_alloc_complex(n);
  Y  = fftw_alloc_complex(n);
  
  // Note that fftw_alloc_complex is equivalent to (fftw_complex *)
  
  x = (float *)malloc(n*sizeof(float));
  f = (float *)malloc(n*sizeof(float));
  
  for(i=0; i < n; i++){
    fscanf(fp1,"%f %f %f\n",&x[i],&t2,&t3);
    X[i][0] = t2;
    X[i][1] = t3;
  }
  fclose(fp1);

  modes(n,f);
  
  dx = (x[n-1]-x[0])/(n-1); 
  f_nyq = 1.0/(2.0*dx);

  // The frequnecy range must be [-f_nyq,f_nyq] which corresponds to [0,n-1] or [-n/2,n/2]
  // note that n/2 corresponds to both +f_nyq and -f_nyq check Numerical Recepies 
  
  for(i=0; i < n; i++)
    f[i] *=2.0*f_nyq;  
  
  p = fftw_plan_dft_1d(n, X,Y, FFTW_FORWARD, FFTW_ESTIMATE); 
  fftw_execute(p); 
  
  q=fftw_plan_dft_1d(n, Y, X1, FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_execute(q);
  
  fp1=fopen("out_1d.dat","w");
  fp2=fopen("rec_1d.dat","w");
 
  for(i=0; i < n; i++){
    fprintf(fp1,"%.6f   %.6e  %.6e\n",f[i],Y[i][0],Y[i][1]);
    fprintf(fp2,"%.6f   %.6e  %.6e\n",x[i],C*X1[i][0],C*X1[i][1]);
  }
  fclose(fp1);
  fclose(fp2);
  

  free(x);
  free(f);
  
  fftw_free(X);
  fftw_free(X1); 
  fftw_free(Y);  
  
  return(0);
  
}