/*  A program to test 3d complex forward and inverse fast fourier transform routines	*/

#include <NR.H>	/* uses fourn from numerical recipes in C to verify ifft3d */
			/*change fmin in numerical recipes to fminnr to avoid conflict with fp.h */

#include <stdio.h>
#include <stdlib.h>
#include <fp.h>
#include <math.h>
#include "fftlib.h"
#include "fftext.h"
#include "fft2d.h"

#if macintosh
#include <timer.h>
#endif

#define NSIZES 	24		/* the number of different ffts col sizes to test */

#define	BIPRAND(a) (2.0/(RAND_MAX+1.0)*a-1.0)

typedef  struct{
	float Re;
	float Im;
	} Complex;

void main(){
long 	fftSize[NSIZES] = 	/* size of FFTs cols, must be powers of 2 */
		{2, 		4, 		8, 		16, 		32, 		64, 		128, 	256,
		512, 	1024, 	2048, 	4096, 	8192, 	16384, 	32768, 	65536,
		131072, 	262144, 	524288, 	1048576, 	2097152, 	4194304, 	8388608, 	16777216};
Complex	*a1;
long 	isize;
long 	i1;
long 	TheErr;
long		N;
long		M;
long		N2 = 16;	/* the number of rows in the 3d fft */
long		M2;
long		N3 = 32;	/* the number of pages in the 3d fft */
long		M3;
float 	maxerrifft;
float 	maxerrfft;
unsigned long 	nn[3];

unsigned int	randseed = 777;
int		rannum;
#if macintosh
	UnsignedWide 		TheTime1;
	Microseconds(&TheTime1);
	randseed = TheTime1.lo;
#endif

printf(" %6d  Byte Floats \n", sizeof(a1[0].Re));
printf(" randseed = %10u\n", randseed);
for (isize = 0; isize < NSIZES; isize++){

	srand(randseed);
	N = fftSize[isize];
	M = roundtol(LOG2(N));
	N = POW2(M);
	M2 = roundtol(LOG2(N2));
	N2 = POW2(M2);
	M3 = roundtol(LOG2(N3));
	N3 = POW2(M3);

	printf("ffts size = %5d X%5d X%6d,  ", N3, N2, N);

	nn[0] = N3;
	nn[1] = N2;
	nn[2] = N;
	
	TheErr = fft3dInit(M3, M2, M);

	if(!TheErr){
		a1 = (Complex *) malloc(N3*N2*N*sizeof(Complex) );
		if (a1 == 0) TheErr = 2;
	}

	if(!TheErr){

			/*  set up a simple test case */
		for (i1=0; i1<N3*N2*N; i1++){
			rannum = rand();
			a1[i1].Re = BIPRAND(rannum);
			rannum = rand();
			a1[i1].Im = BIPRAND(rannum);
		}

			/*  first use fourn from numerical recipes in C to verify ifft3d */
			/*  Note their inverse fft is really the conventional forward fft */
		fourn((float *)a1-1, nn-1, 3, -1);

		ifft3d((float *)a1, M3, M2, M);

		maxerrifft = 0;
		srand(randseed);
		for (i1=0; i1<N3*N2*N; i1++){
			rannum = rand();
			maxerrifft = fmax(maxerrifft, fabs(BIPRAND(rannum)-a1[i1].Re));
			a1[i1].Re = BIPRAND(rannum);
			rannum = rand();
			maxerrifft = fmax(maxerrifft, fabs(BIPRAND(rannum)-a1[i1].Im));
			a1[i1].Im = BIPRAND(rannum);
		}

		printf("errifft = %4.3e,  ", maxerrifft);

			/*  now use iffts to verify ffts */
		ifft3d((float *)a1, M3, M2, M);
		fft3d((float *)a1, M3, M2, M);

		maxerrfft = 0;
		srand(randseed);
		for (i1=0; i1<N3*N2*N; i1++){
			rannum = rand();
			maxerrfft = fmax(maxerrfft, fabs(BIPRAND(rannum)-a1[i1].Re));
			rannum = rand();
			maxerrfft = fmax(maxerrfft, fabs(BIPRAND(rannum)-a1[i1].Im));
		}

		printf("errfft = %4.3e\n", maxerrfft);

		free(a1);
		fft3dFree();
	}
	else{
		if(TheErr==2)	printf(" out of memory \n");
		else	printf(" error \n");
		fft3dFree();
	}
}
printf(" Done. \n");
return;
}