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

int modes(int ,float []);

int main(int argc, char *argv[]) {
  fftw_complex *X, *Y,*X1; 
  fftw_plan p,q;
  int i, j, n;
  FILE *fp1,*fp2,*fp3;
  float *x,*y,*fx,*fy,C1,t1,t2,dx,f_nyq;

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

  fp1=fopen(argv[1],"r");
  n = atoi(argv[2]);
  
  x=(float *)malloc(n*n*sizeof(float));
  y=(float *)malloc(n*n*sizeof(float)); 
  fx =(float *)malloc(n*sizeof(float)); 
  fy =(float *)malloc(n*sizeof(float)); 


  X  = fftw_alloc_complex(n*n);
  X1 = fftw_alloc_complex(n*n);
  Y  = fftw_alloc_complex(n*n);

  for (i=0; i < n; i++){
    for(j=0; j < n; j++){
      fscanf(fp1,"%f %f %f %f\n",&x[j+n*i],&y[j+n*i],&t1,&t2);
      X[j+n*i][0]=t1;
      X[j+n*i][1]=t2;
    }
  }
  fclose(fp1);

  modes(n ,fx);
  modes(n ,fy);

  dx = (x[n*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++){
    fx[i] *=2.0*f_nyq;
    fy[i] *=2.0*f_nyq;
 }
 
  p=fftw_plan_dft_2d(n,n,X,Y, FFTW_FORWARD, FFTW_ESTIMATE); 
  fftw_execute(p); 
  
  q=fftw_plan_dft_2d(n,n,Y,X1, FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_execute(q);
  
  fftw_destroy_plan(p);
  fftw_destroy_plan(q);

  fp1=fopen("out_2d.dat","w");
  fp2=fopen("rec_2d.dat","w");
  
  C1=1.0/(n*n);

  for(i=0; i < n; i++){
    for(j=0; j < n; j++){
      fprintf(fp1,"%.6f %.6f %.6f %.6f\n",fx[i],fx[j],Y[j+n*i][0],Y[j+n*i][1]);
      fprintf(fp2,"%.6f %.6f %.6f %.6f\n",x[j+n*i],y[j+n*i],C1*X1[j+n*i][0],C1*X1[j+n*i][1]);
    }
  }
  fclose(fp1);
  fclose(fp2);

  free(fx);
  free(fy);
  free(x);
  free(y);

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

  return(0);
}