#include <stdlib.h>
#include <math.h>
#include "geostat.h"
#include "toolsludo.h"


/*FAST FOURIER TRANSFORM Pressure Simulation   */
/*Turns a Gaussian white noise vector into a   */
/*spatially correlated vector                  */
/*input:                                       */
/*variogram: structure defining the variogram  */
/*           model                             */
/*grid: structure defining the grid            */
/*n: vector with the number of cells along the */
/*   X, Y and Z axes for the underlying grid   */
/*   i = [0 1 2]                               */
/*   --> 0 0 0        : n will be computed and */
/*   updated as output                         */
/*   --> nx ny nz: these dimensions are used   */
/*realin: structure defining a realization -   */
/*        must be a Gaussian white noise       */
/*gradient: macroscopic gradient pression vector */
/*output:                                      */
/*realout: structure defining a realization -  */
/*realout2: structure defining a pressure field */
/*realout3: structure defining a xvelocity field */
/*realout4: structure defining a yvelocity field */
/*realout5: structure defining a zvelocity field */


void FFTPSim(struct vario_mod variogram,struct statistic_mod stat,struct grid_mod grid,int n[3],struct realization_mod *realin,struct pressure_mod gradient,struct realization_mod *realout,struct realization_mod *realout2,struct realization_mod *realout3,struct realization_mod *realout4,struct realization_mod *realout5)

{
	int NTOT,i,j,k,NMAX,NDIM,ntot,nmax,NXYZ,nxyz,maille0,maille1;
	double *workr,*worki,temp,temp2,coeff;
	double *covar,*realization,*pressure;
	double *icovar,*ireal,*ipressure;
	double *xvelocity,*ixvelocity,*yvelocity,*iyvelocity,*zvelocity,*izvelocity;
	double ki,kj,kk;

	/*covariance axis normalization*/
	axes(variogram.ap,variogram.scf,variogram.Nvario);

	/*pseudo-grid definition*/
	cgrid(variogram,grid,n);  
	
	/*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;

	/*array initialization*/
	covar = (double *) malloc(ntot * sizeof(double));
	testmemory(covar);

	icovar = (double *) malloc(ntot * sizeof(double));
	testmemory(icovar);

	realization = (double *) malloc(ntot * sizeof(double));
	testmemory(realization);

	ireal = (double *) malloc(ntot * sizeof(double));
	testmemory(ireal);

	pressure = (double *) malloc(ntot * sizeof(double));
	testmemory(pressure);

	ipressure = (double *) malloc(ntot * sizeof(double));
	testmemory(ipressure);

	xvelocity = (double *) malloc(ntot * sizeof(double));
	testmemory(xvelocity);

	ixvelocity = (double *) malloc(ntot * sizeof(double));
	testmemory(ixvelocity);

	yvelocity = (double *) malloc(ntot * sizeof(double));
	testmemory(yvelocity);

	iyvelocity = (double *) malloc(ntot * sizeof(double));
	testmemory(iyvelocity);

	zvelocity = (double *) malloc(ntot * sizeof(double));
	testmemory(zvelocity);

	izvelocity = (double *) malloc(ntot * sizeof(double));
	testmemory(izvelocity);

	workr = (double *) malloc(nmax * sizeof(double));
	testmemory(workr);

	worki = (double *) malloc(nmax * sizeof(double));
	testmemory(worki);

	/*covariance function creation*/
	covariance(covar,variogram,grid,n);

	/*power spectrum*/
	fourt(covar,icovar,n,NDIM,1,0,workr,worki);

	/*organization of the input Gaussian white noise*/

	prebuild_gwn(grid,n,realin,realization);

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

	/* build realization in spectral domain   */
	build_real(n,NTOT,covar,realization,ireal);

	free(covar);

      	/* pressure calculation in the spectral domain*/

	build_pressure(n,grid,gradient,realization,ireal,pressure,ipressure);
	build_velocity(n,grid,stat,gradient,realization,ireal,xvelocity,ixvelocity,1);
	build_velocity(n,grid,stat,gradient,realization,ireal,yvelocity,iyvelocity,2);
	build_velocity(n,grid,stat,gradient,realization,ireal,zvelocity,izvelocity,3);

	/*backward fourier transform*/
	fourt(realization,ireal,n,NDIM,0,1,workr,worki);
       	fourt(pressure,ipressure,n,NDIM,0,1,workr,worki);
       	fourt(xvelocity,ixvelocity,n,NDIM,0,1,workr,worki);
       	fourt(yvelocity,iyvelocity,n,NDIM,0,1,workr,worki);
       	fourt(zvelocity,izvelocity,n,NDIM,0,1,workr,worki);
		
	free(ireal);
	free(ipressure);
	free(ixvelocity);
	free(iyvelocity);
	free(izvelocity);
	free(workr);
	free(worki);
	
	/*output realization*/
 	/*is the output realization already allocated?*/
	/*if not, memory allocation*/

	clean_real(realin,n,grid,realization,realout);
	clean_real(realin,n,grid,pressure,realout2);
	clean_real(realin,n,grid,xvelocity,realout3);
	clean_real(realin,n,grid,yvelocity,realout4);
	clean_real(realin,n,grid,zvelocity,realout5);

	free(realization);
	free(pressure);
	free(xvelocity);
	free(yvelocity);
	free(zvelocity);

	return;
}