/**************************************************************************
 *                                                                        *
 * Copyright (c) 1996 Michael Richmond and Richard Treffers               *
 *                                                                        *
 *                    This software may be copied and distributed for     *
 *                    educational, research and not for profit services   *
 *                    provided that this copyright and statement are      *
 *                    included in all such copies.                        *
 *                                                                        *
 **************************************************************************/


/*	
	Numerical recipes version of 2-d transform
	Modified:
		1. for HUGEFLOAT data type
		removed INDEX 9/12/91 RRT
		
	10/29/92 HUGEFLOAT removed
	4/22/94 - ournr.h used instead of nr.h
*/

#include <math.h>
#include "ournr.h"
#include "pcvista.h"

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void 
fourn(data, nn, ndim, isign)
float *data;
int *nn, ndim, isign;
{
	int i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
	int ibit, idim, k1, k2, n, nprev, nrem, ntot;
	float tempi, tempr;
	double theta, wi, wpi, wpr, wr, wtemp;

	ntot = 1;
	for (idim = 1; idim <= ndim; idim++)
		ntot *= nn[idim];
	nprev = 1;
	for (idim = ndim; idim >= 1; idim--) {
		n = nn[idim];
		nrem = ntot/(n*nprev);
		ip1 = nprev << 1;
		ip2 = ip1*n;
		ip3 = ip2*nrem;
		i2rev = 1;
		for (i2 = 1; i2 <= ip2; i2 += ip1) {
			if (i2 < i2rev) {
				for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2) {
					for (i3 = i1; i3 <= ip3; i3 += ip2) {
						i3rev = i2rev +i3 - i2;
						SWAP(data[i3], data[i3rev]);
						SWAP(data[i3 + 1], data[i3rev + 1]);
					}
				}
			}
			ibit = ip2 >> 1;
			while (ibit >= ip1 && i2rev > ibit) {
				i2rev -= ibit;
				ibit >>= 1;
			}
			i2rev += ibit;
		}
		ifp1 = ip1;
		while (ifp1 < ip2) {
			ifp2 = ifp1 << 1;
			theta = isign*6.28318530717959/(ifp2/ip1);
			wtemp = sin(0.5*theta);
			wpr = -2.0*wtemp*wtemp;
			wpi = sin(theta);
			wr = 1.0;
			wi = 0.0;
			for (i3 = 1; i3 <= ifp1; i3 += ip1) {
				for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2) {
					for (i2 = i1; i2 <= ip3; i2 += ifp2) {
						k1 = i2;
						k2 = k1 + ifp1;
						tempr = (float) (wr*data[k2] - wi*data[k2 + 1]);
						tempi = (float) (wr*data[k2 + 1] + wi*data[k2]);
						data[k2] = data[k1] - tempr;
						data[k2 + 1] = data[k1 + 1] - tempi;
						data[k1] += tempr;
						data[k1 + 1] += tempi;
					}
				}
				wr = (wtemp = wr)*wpr - wi*wpi + wr;
				wi= wi*wpr + wtemp*wpi + wi;
			}
			ifp1 = ifp2;
		}
		nprev *= n;
	}
}

#undef SWAP