/*	Ft2d.c

	2-d Fourier transform subroutine using the 1-d FT subroutine.

	Usage:		int n1, n2;
			float sign1, sign2;
			struct complex cp[][], cwork[];
			ft2d(n1, n2, cp, sign1, sign2, cwork);

	Compile:	cc prog.c ft2d.c fft.c -fswitch -lm -o prog

	Arguments:	n1	number of rows of 2-d matrix
				MUST BE A POWER OF 2
			n2	number of columns of matrix, fastest varying
				MUST BE A POWER OF 2
			cp	pointer to complex 2-d matrix cp[n1][n2]
			sign1	sign for row transform
			sign2	sign for column transform
			cwork	pointer to scratch vector of n1 complex
				elements
									*/
/* Include the usual header files */
#include <stdio.h>
#include <math.h>

/* Define complex multiplication and its conjugate */
#define  rmul(x,y)      (x->re * y->re - x->im * y->im)
#define  imul(x,y)      (x->im * y->re + x->re * y->im)
#define rcmul(x,y)      (x->re * y->re + x->im * y->im)
#define icmul(x,y)      (x->im * y->re - x->re * y->im)

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

/* Body of subroutine */
ft2d(n1, n2, cp, sign1, sign2, cwork)
/* Declare all argument variables */
int n1, n2;
float sign1, sign2;
struct complex *cp, *cwork;
	{
	/* Declare all local variables */
	int i1, i2;
	FILE *rep;
	rep = fopen("report", "a");

	/* Transform over the slow dimension, the columns */
	for(i2=0; i2<n2; i2++) /* Loop through columns */
		{
		/* Copy each column into scratch vector */
		for(i1=0; i1<n1; i1++) /* Loop through rows */
			{
			cwork[i1].re = cp[i1*n2+i2].re;
			cwork[i1].im = cp[i1*n2+i2].im;
			}
		/* One-dimensional Fourier transform */
		fft(n1, cwork, sign2, 1.0);
		/* Copy the columns back into the matrix */
		for(i1=0; i1<n1; i1++) /* Loop through rows */
			{
			cp[i1*n2+i2].re = cwork[i1].re;
			cp[i1*n2+i2].im = cwork[i1].im;
			}
		if (i2==0 || i2==1 || i2%20 == 0)
			{
			fprintf(rep, "transformed column %d\n", i2);
			fflush(rep);
			}
		}

	/* Transform over the fast dimension, the rows */
	for(i1=0; i1<n1; i1++)
		{
		fft(n2, &(cp[i1*n2+0]), sign1, 1.0);
		if (i1==0 || i1==1 || i1%20 == 0)
			{
			fprintf(rep, "transformed row %d\n", i1);
			fflush(rep);
			}
		}

	fclose(rep);
	/* Return to main program */
	return(0);
	}