/*	Fttest.c

	Test case for two-dimensional Fourier transform.
	From Claerbout (1985) p. 69.

	Usage:		fttest >outfile

	Compile:	cc fttest.c ft2d.c fft.c -lm -o fttest
								*/

/* Get the usual standard output and math definitions */
#include <stdio.h>
#include <math.h>

/* Declare the structure to hold complex numbers */
struct complex 
	{
	float re; 
	float im;
	};

/* Enter main part of program */
main(ac, av)
int ac;
char **av;
	{
	/* Declare all variables */
	int it, nt, ix, nx;
	/*
	float p[64][64];
	struct complex cp[64][64], cwork[64]; * Use rows of time *
	*/
	float *p, *cwork;
	struct complex *cp;

	/* Assign sizes */
	if (ac < 2 || ac > 3)
		nx = 512;
	else
		nx = atoi(av[1]);
	if (ac == 3)
		nt = atoi(av[2]);
	else
		nt = nx;

	cp = (struct complex *)malloc(nx*nt*sizeof(struct complex));
	allocerr(cp, nx*nt*sizeof(struct complex), "cp");
	p = (float *)malloc(nx*nt*sizeof(float));
	allocerr(p, nx*nt*sizeof(float), "p");
	cwork = (float *)malloc(nx*nt*sizeof(float));
	allocerr(cwork, nx*nt*sizeof(float), "cwork");

	/* Nested loops to set complex array to zero; time varies fastest */
	for(ix=0; ix<nx; ix++)
		{
		for(it=0; it<nt; it++)
			{
			cp[ix*nt+it].re = 0;
			cp[ix*nt+it].im = 0;
			}
		}

	/* Insert a spike into the real part of the complex array */
	cp[3*nt+16].re = 1; cp[4*nt+16].re = 4; cp[5*nt+16].re = 6; cp[6*nt+16].re = 4;
	cp[7*nt+16].re = 1; cp[3*nt+17].re = 1; cp[4*nt+17].re = 4; cp[5*nt+17].re = 6;
	cp[6*nt+17].re = 4; cp[7*nt+17].re = 1;

	/* Transform array into the Fourier domain */
	ft2d(nx, nt, cp, 1.0, 1.0, cwork);

	/* Copy the real part of the transformed complex array to a real array
	   to allow writing */
	for(ix=0; ix<nx; ix++)
		{
		for(it=0; it<nt; it++)
			{
			p[ix*nt+it] = cp[ix*nt+it].re;
			}
		}
	/* Write 2-d real array of transform in binary to standard output */
	write(1, p, nx*nt*sizeof(float));

	/* All done */
	exit(0);
	}

#define NAMEL 120
#define EROUT stderr

allocerr(p, nbytes, s)			/* error handling for malloc */
int p, nbytes;
char s[NAMEL];
	{
	if(p==NULL)
		{
		fprintf(EROUT, "allocerr: cannot allocate %d bytes of memory for %s\n",
		nbytes, s);
		if (nbytes > 100000000)
			fprintf(EROUT, "Don't be so greedy!\n");
		exit(-1);
		}
	return(0);
	}