//This example code is used
//to familiarize one with what exactly
//and how FFTW performs FFT.
// compile by
// g++ -lm -I/usr/local/include -L/usr/local/lib -lrfftw -lfftw fftw_exercise_3by3.C
// The routine is from fftw version 2. The newest version 3 might have different
// format.

#include <fstream>
#include <iostream>
#include "rfftw.h"

using namespace std;

const int  Nmesh=3;
const int  Nmesh2=9;
const double twopi=6.283185307;
 
int main(int argc, char **argv)
{
  char   buf[200];
  int    nyquist[3],blklen;
  int    ratio, axes, zone,i,j;
  double kmin, kmax;
  float xdummy;
  float kvector[3];
  char inname[100],outname[100];

  fftw_complex *A, *B;
  rfftwnd_plan p, pinv;
  fftw_real scale = 1.0 / (Nmesh *Nmesh);
 

  if(argc<1)
    {
      cout<<"Parameters are missing"<<endl;
      cout<<"Call with <ParameterFile>"<<endl;
    }

  fftw_real * real_array = new fftw_real[Nmesh*(2*(Nmesh/2+1))];

  //the input array is, in row-major order

  //   A(0,0)  A(0,1)  A(0,2)  =   1   2   3
  //   A(1,0)  A(1,1)  A(1,2)  =   4   5   6
  //   A(2,0)  A(2,1)  A(2,2)  =   7   8   9
  
  //initialise
  //  for(i=0;i<Nmesh;i++)
  //   for(j=0;j<Nmesh;j++)
  //      real_array[i * (2*(Nmesh/2+1)) + j]=xdummy; //padding fftw style

  
  i=0; j=0;
  real_array[i * (2*(Nmesh/2+1)) + j]=1;

  i=0; j=1;
  real_array[i * (2*(Nmesh/2+1)) + j]=2;

  i=0; j=2;
  real_array[i * (2*(Nmesh/2+1)) + j]=3;

  i=1; j=0;
  real_array[i * (2*(Nmesh/2+1)) + j]=4;

  i=1; j=1;
  real_array[i * (2*(Nmesh/2+1)) + j]=5;

  i=1; j=2;
  real_array[i * (2*(Nmesh/2+1)) + j]=6;

  i=2; j=0;
  real_array[i * (2*(Nmesh/2+1)) + j]=7;

  i=2; j=1;
  real_array[i * (2*(Nmesh/2+1)) + j]=8;

  i=2; j=2;
  real_array[i * (2*(Nmesh/2+1)) + j]=9;


  p    = rfftw2d_create_plan(Nmesh, Nmesh, FFTW_REAL_TO_COMPLEX,
                             FFTW_ESTIMATE | FFTW_IN_PLACE);
  pinv = rfftw2d_create_plan(Nmesh, Nmesh, FFTW_COMPLEX_TO_REAL,
                             FFTW_ESTIMATE | FFTW_IN_PLACE);
  
     
  cout<<"FFT..."<<endl;
  rfftwnd_one_real_to_complex(p, &real_array[0], NULL);
  
  // aliases for accessing complex transform outputs: 
  A = (fftw_complex*) &real_array[0];

  for(i=0;i<Nmesh;i++){
    for(j=0;j<Nmesh/2+1;j++){
    //for(j=0;j<Nmesh;j++){
      cout<<i<<" "<<j<<" "<<" ("<<i*(Nmesh/2+1)+j<<") "<<A[i*(Nmesh/2+1)+j].re<<" "<<A[i*(Nmesh/2+1)+j].im<<endl;
    }
  }
  
  cout<<"Inverse FFT..."<<endl;
  rfftwnd_one_complex_to_real(pinv, A, NULL);
  

  for(i=0;i<Nmesh;i++){
    for(j=0;j<Nmesh;j++){
        
      cout<<i<<" "<<j<<" "<<real_array[i * (2*(Nmesh/2+1)) + j]/(double)(Nmesh*Nmesh)<<endl;

    }
  }

  
  rfftwnd_destroy_plan(p);
  rfftwnd_destroy_plan(pinv);
  
}