#include<stdio.h>
#include<stdlib.h>
#include<fftw3.h>
#include<math.h>
#include<complex.h>
#include <omp.h>

int modes(int ,float []);

int main (int argc, char *argv[]){
  fftw_complex *X, *X1, *Y;  
  fftw_plan p,  q; 
  int i, j, k,l,i1,j1,k1, n;
  FILE *fp1,*fp2;
  float *x,*y,*z,*fx,*fy,*fz,C,t1,t2,f_nyq,dx; 

  if (argc < 2){
    fprintf(stderr,"./example_dft_3d  <input file> <size along one dim>\n");
    return(-1);
  }

  fp1=fopen(argv[1],"r");
  n=atoi(argv[2]);

  printf("max threads=%d\n",omp_get_max_threads());

  fftw_init_threads();
  fftw_plan_with_nthreads(omp_get_max_threads());
  //fftw_plan_with_nthreads(8);

  X  = fftw_alloc_complex(n*n*n);
  X1 = fftw_alloc_complex(n*n*n);
  Y  = fftw_alloc_complex(n*n*n);
 
  x=(float *)malloc(n*n*n*sizeof(float));
  y=(float *)malloc(n*n*n*sizeof(float));
  z=(float *)malloc(n*n*n*sizeof(float));

  fx=(float *)malloc(n*sizeof(float));
  fy=(float *)malloc(n*sizeof(float));
  fz=(float *)malloc(n*sizeof(float));

  for(i=0; i < n; i++){
    for(j=0; j < n; j++){
      for(k=0; k < n; k++){
        l = k+n*(j+n*i); 
	fscanf(fp1,"%f %f %f %f %f\n",&x[l],&y[l],&z[l],&t1,&t2);
        X[l][0] = t1; 
        X[l][1] = t2; 
      }
    }
  }
  fclose(fp1);
 
  modes(n ,fx);
  modes(n ,fy);
  modes(n ,fz);
  C = 1.0/(n*n*n); 

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

 printf("f_nyq=%.6f\n",f_nyq);

  // 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++){
    fx[i] *=2.0*f_nyq;
    fy[i] *=2.0*f_nyq;
    fz[i] *=2.0*f_nyq;
 }

  p=fftw_plan_dft_3d(n,n,n, X, Y, FFTW_FORWARD, FFTW_ESTIMATE);  
  fftw_execute(p); 
  
  // This is the inverse fourier transform 
    
  q=fftw_plan_dft_3d(n,n,n, Y, X1, FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_execute(q); 
  
  fp1=fopen("out_3d.dat","w");
  fp2=fopen("rec_3d.dat","w");
  
  for(i=0; i < n; i++){
    for(j=0; j < n; j++){
      for(k=0; k < n; k++){
        l = k+n*(j+n*i);   
	fprintf(fp1,"%.6f %.6f %.6f %.6e %.6e\n",fx[i],fy[j],fz[k],Y[l][0],Y[l][1]);
	fprintf(fp2,"%.6f %.6f %.6f %.6f %.6f\n",x[l],y[l],z[l],C*X1[l][0],C*X1[l][1]);
      }
    }
  }
 
  free(x);
  free(y);
  free(z);

  free(fx);
  free(fy);
  free(fz); 

  fftw_free(X);
  fftw_free(Y);
  fftw_free(X1);



  fclose(fp1);
  fclose(fp2);
  
  return(0);
  
  
}