#include <stdio.h>
#include <stddef.h>
#include <string.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include "geostat.h"


/*FAST FOURIER TRANSFORM Test   */

void FFTtest(int n[3],struct grid_mod grid,struct realization_mod *realin,struct realization_mod *realout,int solver)

{
  int NTOT,i,j,k,NMAX,NDIM,ntot,nmax,NXYZ,nxyz,maille0,maille1;
  double *workr,*worki,temp,temp2,coeff;
  double *realization;
  double *ireal;
  double ki,kj,kk;
  FILE *fp;

  /*constant definition*/
  NTOT = n[0]*n[1]*n[2];
  ntot = NTOT+1;
  NMAX = n[0];
  NDIM = 3;
  for (i=1;i<3;i++) {
    if (n[i] > NMAX)
      NMAX = n[i];
    if (n[i] == 1)
      NDIM--;
  }
  nmax = NMAX+1;
  NXYZ = grid.NX*grid.NY*grid.NZ;
  nxyz = NXYZ+1;
  
  realization = (double *) malloc(ntot * sizeof(double));
  testmemory(realization);
  
  ireal = (double *) malloc(ntot * sizeof(double));
  testmemory(ireal);
  
  workr = (double *) malloc(nmax * sizeof(double));
  testmemory(workr);
  
  worki = (double *) malloc(nmax * sizeof(double));
  testmemory(worki);
  
  /*organization of the realization*/
  prebuild_gwn(grid,n,realin,realization,solver);

  fp = fopen("perm.test", "w");
  for (i=1;i<=NTOT;i++) 
    fprintf(fp,"%f\n",realization[i]);
  fclose(fp);

  /*forward fourier transform of the GWN*/
  fourt(realization,ireal,n,NDIM,1,0,workr,worki);
  
  /*backward fourier transform*/
  fourt(realization,ireal,n,NDIM,0,1,workr,worki);

  fp = fopen("perm2.test", "w");
  for (i=1;i<=NTOT;i++) 
    fprintf(fp,"%f\n",realization[i]);
  fclose(fp);
  return;
}