/*	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 -fswitch -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()
	{
	/* Declare all variables */
	int it, nt, ix, nx;
	float p[64][64];
	struct complex cp[64][64], cwork[64]; /* Use rows of time */

	/* Assign sizes */
	nx = 64;
	nt = 64;

	/* 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][it].re = 0;
			cp[ix][it].im = 0;
			}
		}

	/* Insert a spike into the real part of the complex array */
	cp[3][16].re = 1; cp[4][16].re = 4; cp[5][16].re = 6; cp[6][16].re = 4;
	cp[7][16].re = 1; cp[3][17].re = 1; cp[4][17].re = 4; cp[5][17].re = 6;
	cp[6][17].re = 4; cp[7][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][it] = cp[ix][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);
	}